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I . INTRODUCTION 

Under  previous  contracts  (N61331-74-C-0027  and  N61331-74-C-0070) 
and  the  present  contract  (N61331-75-C-0012)  a new  technique  for  the 
Identification  of  vehicle  transfer  functions  from  flight  test  data  has 
been  developed  and  tested.  The  Importance  of  dynamics  Identification 
arises  because  the  resulting  model  provides  1)  a quantitative  measure 
of  the  stability  of  the  vehicle,  11)  a means  to  predict  the  vehicle 
response  to  new  maneuvers,  and  111)  a basis  for  designing  and  evaluating 
a control  system  for  the  vehicle,  through  simulation  or  analytical  pro- 
cedures. To  elaborate  on  the  latter  consider  a towed  sonar  vehicle  for 
which  the  deviations  of  the  roll,  pitch  and  yaw  angles  — l.e.  the  vehicle 
attitude  angles  — are  required  to  be  controlled  carefully.  Clearly,  the 
dynamics  of  the  vehicle  must  first  be  Identified  before  the  analytical 
design  of  an  optimum  controller,  or  the  subsequent  evaluation  of  the  over- 
all system  through  simulation,  can  be  carried  out. 

Salient  features  of  the  new  Identification  technique  GRAM  are 
summarized  below: 

(a)  Noniterative  - The  method  is  non-lteratlve  [1],  and  therefore 
avoids  the  problems  of  convergence  often  associated  with  Iterative 
techniques . 

(b)  Noiseworthy  - The  method  Is  noiseworthy  [2].  Accurate 
estimation  can  therefore  be  achieved  even  In  the  presence  of 
high  levels  of  noise  [2]  provided  that  the  system  response  con- 
tains all  of  its  modes  In  a reasonably  observable  manner.  This 
success  can  be  attributed  to  the  fact  that  the  method  utilizes 
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signals  generated  via  measurement  filters  whose  poles  can  be 
chosen  appropriately  by  the  test  engineer. 
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automatic,  such  that  the  overall  system  behaves  In  a stable  manner. 
Then  It  Is  found  that  the  application  of  GRAM  technique  to  the  Input- 
output  data  of  the  unstable  vehicle  yields  its  dynamics,  including 
the  poles  In  the  right-half  s-plane,  with  considerable  accuracy. 


Under  the  present  contract,  #N61331-75-C-0012,  the  following  segments  of 
research  work  were  performed. 

1.  Consideration  of  non-zero  Initial  conditions  and  state  model 
formulation. 

2.  Application  of  GRAM  technique  to  flight  test  data  of  the  DSRV  (deep 
submergence  research  vehicle) . A special  feature  of  this  work  was  the  detect- 
ion of  a weak  oscillatory  mode  from  the  short  duration  test  data  that  was 

ava I lable . 

3.  Adaptation  and  testing  of  GRAM  technique  to  unstable  systems. 

4.  Cascading  of  the  GRAM  technique  with  the  NASA  (Taylor)  Systems 
Identification  Technique  [3]. 

5.  Application  of  GRAM  to  flight  test  data  generated  from  the  Navy 
standard  program  'DYNAMIC'  In  the  presence  of  mild  nonlinearities. 
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The  above  research  efforts,  which  encompass  theoretical  developments, 

I 

computer  implementation,  and  testing,  are  described  in  the  succeeding  chapters.  j 
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li.  MULTIVARIABLE  IDENTIFICATION  OF  VEHICLE  DYNAMICS  WITH 
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NON-ZERO  INITIAL  CONDITIONS,  STATE-MODEL 

The  purpose  of  this  chapter  Is  Co  discuss  the  identification  of  vehicle 
dynamics  under  non-zero  initial  conditions,  estimating  these  unknown  initial 
conditions,  and  to  generate  a state-variable  model  from  the  identified  transfer 
function  matrix.  A computer  program  written  in  FORTRAN  IV  was  developed  and 
tested.  A listing  of  the  complete  program,  MGRAM,  and  a card  deck  for  the 
same  are  enclosed.  After  a theoretical  discussion  of  the  method  the  results 
of  id.>ntlfication  performed  on  the  simulated  flight  test  data  of  a six-man 
submersible  are  presented. 

2.1  COMMON  MODE  FORMULATION 

Consider  that  the  vehicle  dynamics  Is  described  by  the  state  equations 

(la) 

v(k)  = F^x(k)  + F^u^(k)  (lb) 


x(k+l)  = F^x(k)  + F^u^(k) 


where  x(k)  is  the  N-dimensional  vector  of  state  variables,  v(k)  is  the 

m-dimensional  vector  of  outputs,  and  u^(k)  is  an  r-dimensional  vector  of 

12  3 4 

known  inputs.  The  constant  matrices  F , F , F and  F have  appropriate 
dimensions.  In  general,  the  measured  outputs  are  corrupted  by  bias  and  noise: 

y(k)  = v(k)  + f + e(k)  (Ic) 

The  z-domain  description  of  this  syste-  is  given  by 

x(z)  = (zI-F^)  ^[F^u^(z)  + x(0)z]  (2a) 

y(z)  = F^x(z)  + F^u^(z)  + f (1-z  ^)  ^ + e(z)  (2b) 

or,  upon  combining  these  equations, 

y(z)  * F^(zT-F^)  ^[F^u^(z)+x(0)z]  + F^u^(z)  + f (1-z  S 

+ e(z)  (1) 

where  I denotes  the  identity  matrix  of  order  N.  Upon  rearraiiKemenl  (3)  becomes 


y(z) 


Ih’(z)  ; h(z)  ; f,i 


m'(z) 

1 

L(i-z'‘)"‘j 


+e(z) 


(4) 
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In  (4) 


I 

y 


h\z) 

= 

F^(zT-F^)'V  + F^ 

(5) 

h(z) 

= 

zF^(zI-F^)~^x(0) 

(b) 

The  important 

observation  arises  that 

the  case  of  non-zero  initial 

conditions  has  been 

cast  into  a form  that  is 

analogous  to  the 

case  of  zero 

initial  conditions 

and 

zero  bias.  We  now  define 

H(z) 

= 

[H^(z)  ; h(z)  ; f) 

(7) 

u(z) 

= 

(u^z)'^  ; 1 ; (i-z) 

-IjT 

(8) 

whire  the  superscript  T 

' denotes  the  transpose  operation.  Of 

course,  if  the 

initial  conditions 

and 

the  bias  are  zero  then  the  formulation 

becomes 

= 

n\z) 

(9) 

ii(z) 

= 

u^(z) 

(10) 

Regardless  of  the  above 

details  we  can  now  state  the  problem 

as  follows. 

Problem  statement  - 

Consider  the  model  (see 

Fig.  1 ) 

y (z) 

H(z)  u(z)  + e(z) 

(11) 

where  H(z)  is  an  mxr 

matrix  with  rational 

entries  having  a 

common  denomi 

nutor;  i.e.  its  entries 

have  the  form 

h . . ( z ) 
ij 

Bjj  (z)  / A(z) 

-N. 

. . . + b . . z 

(12) 

°l.iO  "ijl 

1.1  Nj 

(13) 

...  . -N 

N' 

e(k) 


u(k) 


y(k) 


Fig.  1.  Vehicle  model 
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Given  the  model  (11),  the  sequences  y(k)  and  u(k)  over  a set  of  points  | 

0,1,  . . ,L  determine  the  numerator  parameters  1*1,  ..,m,  j-1,  ..,r, 

n“l,  ..,N  and  the  denominator  parameters  a^,  n*!,  . . ,N  such  that  the 
estimated  model  provides  a best  fit  in  some  sense. 

The  error  vector  e(k)  may  arise  from  measurement  errors,  insufficiently  ^ 

high  order  of  approximation,  or  simply  incorrect  modeling  of  the  process. 

Suitable  assumptions  must  therefore  be  made  on  the  error  vector  e(k)  depending  I 

on  the  specific  application. 


2.2  DEVELOPMENT  OF  A KEY  EQUATION 

Before  dealing  with  the  problem  stated  above  let  us  point  out  that 
this  formulation  emphasizes  common  modes  between  the  various  outputs.  This 
is  in  contrast  to  the  formulation  in  [4]  where  each  output  was  considered  to 
have  its  own  denominator  dynamics.  That  freedom  leads  to  undesirable  results 
when  it  is  known  a priori  that  the  outputs  (i.e.,  the  different  entries  of 
the  output  vector  y(k))  share  common  modes.  Such  is  the  case,  for  example, 
when  y(k)  = (u(k) , 0(k) , w(k) 1 , the  longitudinal  variables  of  a submersible. 

Indeed,  we  have  assumed  H(z)  in  (11)  to  have  the  common-mode  form 
H(2)  . tB^j(z)/A(z)]^^^ 

We  now  proceed  with  the  development  by  first  multiplying  equation  (11) 
by  A(z).  This  yields 

A(z)y(/)  = ll(7.)u(z)  + A(z)e(z)  (IM 
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It  Is  Instructive  to  look  at  equation  (14)  explicitly  for  the  two-input, 
I two-output  case.  For,  then. 


1 

A(z)yj^(z) 

B,i(z) 

1 

A(z)y2(z) 

B22^(z) 

Uj^(z) 


U2(z) 


A(z)ej^(z) 


A(z)e2(z) 


(15) 


which,  in  the  time-domain  can  be  written  as 


For  the  general  multi-input,  multi-output  relation  in  (14)  we  define 


V = E (N  +1) 

j=l  ^ 

A = [ 1 


N+1  column  vector 


= 

[bijO  ••• 

N^^l 

column  vector 

= 

[B(il)T  B<^2)T 

g(lr)TjT 

V 

column  vector 

B 

b(2)T  ...  B 

(m)TjT 

mv 

column  vector 

y<^>(k) 

= [yj^(k)  y^(k-i) 

...  y^(k-N)] 

N+1 

row  vector 

e^^^k) 

= [e^(k)  e^(k-l) 

..  e^(k-N)] 

N+1 

row  vector 

u^^^k) 

= [Uj(k)  Uj(k-l) 

..  Uj(k-Nj)] 

N^+1 

row  vector 

U(k)  = 

[u^^^(k)  u^^^(k) 

...  u^’^^k)] 

V 

row  vector 

Using  the  above  notation  we  can  write  (14) 
in  the  time  domain  as  m equations  of  the  form 
y^^^(k)A  = U(k)B^^^  + e^^^(k)A 
where  i=l,  2 ...  m,  or  as  one  equation  as  follows: 


(17) 


-U(k)  0 ...  0 

fA- 

’e^^^k). 

rM 

y<^><k) 

O 

o 

— 

e^^^(k)l 

1 

1 

\ * 

j 

• • • 

S 

• ! 0 

I 

• • • 

B 

. 1 

B 

_y'">(k) 

0 0 ...  -U(kj 

e^^^k)': 

1 

\ 

i 

— 

L J 

(18) 


Denote  the  m x (N+l+mv) dimensional  matrices  on  the  left  and  right  hand  sidis 
of  equation  (18)  as  S and  E,  respectively,  and  the  (N+l+mv)d Imens ioiia I 

rj.  fj»  rp 

parameter  vector  | A | B ] as  Y • Tlien,  we  can  ri-wriie  ( I H)  as 

S(k)Y  = E(k)  Y (''7) 


8 


Note  that  only  the  first  (N+1)  columns  of  E can  have  nonzero  entries; 
the  remaining  columns  are  zero.  In  fact,  we  can  express  E as 

E(k)  - [ Ej(k)l  0 1 (20) 

where  0 is  an  m x V dimensional  matrix  of  zeros. 

Equation  (19)  is  the  key  equation  we  wished  to  obtain  in  this  section. 

The  development  of  the  method, presented  in  the  succeeding  sections,  is  based 
upon  this  relation. 


2 . 3 SOLUTION 

Difficulty  in  solving  (19)  for  the  parameter  vector  y arises  because 
the  error  vector  e(k) , and  therefore  the  matrix  E(k) , is  generally  unknown. 
Suitable  assumptions  on  e(k)  must  therefore  be  made.  We  shall  solve  (19)  for 
increasingly  less  restrictive  assumptions  on  e(k). 

2.3.1  Ideal  Data  (e(k)®0) 

When  the  error  in  observations  y(k)  can  be  assumed  to  be  zero,  equation 
(19)  becomes 

S(k)  Y “ 0,  k»0,l,  ..,  K (21) 


The  innerproduct  of  (21)  when  summed  over  all  k yields 


where 


y"  G Y = 0 

K „ 

G = T.  S (k)  S(k) 
k=0 

is  the  Gram  matrix  of  measurement  signals.  Note  that  (!  Is  real  symmul r ic , 


(22) 

(21) 


hence  all  of  its  eigenvalues  are  real  and  nonnegative.  Equation  (22)  indicates 

that  one  of  the  eigenvalues  is  zero;  call  the  corresponding  eigenvector  as 

..(0) 


Then  the  desired  parameter  vector  is  given  by 

/O) 

Y m I 


(24) 
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2.3.2  Additive  White  Noise  In  Observed  Signals 


Suppose  the  error  in  observations  can  be  assumed  to  be  zero  mean 
white  noise,  i.e. , 


E [e(k)I  - 0 

I [e(k)  eW]  = D 


where  E denotes  the  expectation  operator,  D a diagonal  matrix,  and  6 
the  usual  Knoecker  delta.  Let  the  diagonal  entries  of  D be  denoted  as  ‘‘ii 
and  the  trace  as  d.  Then  the  inner  product  of  (19)  with  itself  when 
summer  over  k=0,  ...,  K yields 


Here, 


Note  that 


T T 

Y G Y = Y P Y 


K 


E P 


z S"(k)  S(k), 
k=0 

T 

Z E (k)  E(k) 
k=0 


K m 

I E R 

k=0  1-1 

-^(11) 


= BC  = 6 


C’ 

K) 


(i) 

0 

0 


(25) 

(26) 

(27) 

(28) 

(29) 


where  is  the  correlation  matrix  of  the  ith  row  of  the  matrix 

E(k) , 6 = (K+l)d  and  is  an  identity  matrix  of  order  N+1.  We  shall 

solve  the  problem  of  estimating  y as  that  of  solving 


Y Gy  - By  p cy  no) 

such  that  B is  minimized.  (Note  that  for  3=0  the  problem  reduces  to 

case  A.)  That  the  estimate  so  obtained  makes  sense  is  justifiable  in  two  ways. 

• When  the  number  of  observations  is  large,  P is  expected  to  tend 

to  3C;  thu.s,  in  the  absence  of  the  knowledge  of  P we  may  replate 
it  by  B<;  In  e(|iiatlon  (25)  to  obtain  (30). 

• If,  instead  of  using  successive  values  of  k in  the  summations 


10 


(26)  and  (27),  we  use  k»0,  N+1,2(N+1),  ...  then  it  can  be  shown 
that  (30)  provides  a maximum  likelihood  estimate  [16]  lor  the 
parameter  vector  under  the  assumption  of  gausslan  white  noise. 

Since  a non-zero  solution  Is  required,  solving  (30)  is  clearly 
equivalent  to  solving  the  generalized  eigenvalue  problem 


Gy  = 6 C Y 
or,  in  partitioned  form. 


for  the  minimum  eigenvalues  and  the  corresponding  eigenvector.  The  above 
problem  is  equivalent  to  solving  the  pair 


[(Gil  - ^12^22^^21^  ~ ^ ^ ^ 


,(2)  , _ g -Ig  Y 

r » - ^22  ^2i 


The  first  part  is  solved  as  a usual  eigenvalue  problem.  The  eigenvector 

Y corresponding  to  the  minimum  eigenvalue  is  selected  and,  then,  from 

(2) 


the  second  equation  Y 
finally  obtained  as 


is  obtained.  The  desired  parameter  vector  is 


"I 


.(2) 


(31) 


(32) 


( } Ja) 


(33b) 


(3M 


2.3.3  Additive  Co  lor ed  Noise  In  Output  Signals 

Suppose  tile  error  in  output  observations  is  assumed  to  be  .i  ,s( .ii  i oii.irv 

random  process  as  shown  in  Figure  2,  where  /a  is  an  unknown  scalar,  l(z)  a 

known  (non-antlclpatory)  transfer  function  matrix  and  w(k)  a white  noise 

T 

vector  sequence;  that  is  £ w(k)w  (2.)  ■ 


11 


z'  {F(2)}-  f(k) 


Fig.  2.  Error  Process 


Using  f(k),  the  matrix  of  impulse  responses  of  the  entries  of  F(z), 
we  can  write 

OC 

e(k)  = /a  Z iil)  w(k-2,) 

!l=0 

And,  it  can  be  shown  easily  that 

a 

E e(k)e^(k+j)  = a E f (il)f^()2.+j)  = OiR(j)  (by  definition) 
il=0 


(3b) 


In  the  interest  of  avoiding  degeneracy  we  assume  that  R(0)  is  positive  definite. 

We  now  return  to  the  main  problem  of  estimating  Y-  For  convenience,  lei 
us  rewrite  equation  (19), 

S(k)Y  = K(k)Y  ^ 

Recall  that  only  the  first  N+1  columns  of  E(k)  are  nonzero.  In  fact  these- 

columns  are  e(k) e(k-N)  respectively.  As  in  subsection  2.3.2,  we  form 

the  outer-product  of  (19)  with  Itself  and  sum  over  k=0,  ..,K.  This  yields  an 
equation  which  is  similar  to  (25)  where,  again,  G and  P are  matrices  given  by 
(26)  and  (27)  respectively.  The  difference  arises  in  the  fact  that  now 
E P = a C 

,(11) 


a 


0 


(3 


0 


0 


where  tlie  (N+l)  x (N+1)  matrix  C 


.(11) 


is  given  .'is 


I 


H 

0 

i: 


i 


,.(11) 


i.i 


m 

>’.  R 
J.=  l 


9J, 


(18) 


Notf  lli.-U  this  mntrJx  Is  symmetric  Toeplitz  and  nonsingiilar . 

Following  a procedure  similar  to  that  in  subsection  4.2,  it  can  be  sitown 
that  the  parameter  vector  y can  be  estimated  by  solving  the  pair 


11 


,(2) 


-G22  Y 


(1) 


(39a) 

(39b) 


The  first  part  is  solved  as  a usual  eigenvalue  problem.  The  eigenvector 

corresponding  to  the  minimum  eigenvalue  is  selected  and,  then,  from 
(2) 

the  second  equation  Y Is  obtained.  The  desired  parameter  vector 
I'stiniate  Is  finally  obtained  as  in  equation  (34)  by  adjoining  these  two  parts 
and  dividing  by  the  first  entry  of  the  first  part. 


4.4  Additive  White  Noise  in  Both  Input  and  Output  Measurements 

The  measured  input,  u(k),  is  considered  to  be  composed  of  measurement 
noise  v(k)  and  the  true  input  u(k)-v(k)  see  Fig.  3.  The  process  v(k)  is 
assumed  to  be  zero  mean  and  white.  The  model  described  in  (11)  must  now 
be  modified  to 

y(z)  = H(z)[u(z)  - v(z)]  + e(z)  (40) 

wlieri-  the  other  symbols  have  their  earlier  defined  meanings.  Upon, 
multiplying  this  equation  by  A(z) , the  common  denominator  of  the  entries 


1 


U 


(41) 


Fig.  3.  Noisy  input  and  output  measurements 
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of  H(z),  one  obtains 


A(z)y(z)  - B(z)u(z)  - B(z)v(z)  + A(z)e(z) 


We  now  claim  that  (36)  can  be  written  in  the  time  domain  as 


S(k)Y  - E(k)Y 

where  S(k)  is  formed  using  y(k),  ...  ,y(k-N)  and  u(k),  ...  ,y(k-N) 
exactly  as  described  in  Section  III.  The  matrix  E(k)  is,  however,  formed 
differently  than  the  E(k)  of  equation  (19).  In  fact,  E(k)  is  now  formed 
by  replacing  y(k)  by  e(k)  and  u(k)  by  v(k)  in  the  matrix  S(k).  Thus 
whereas  E(k)  of  (19)  had  a structure  very  different  from  that  of  S(k), 
the  matrix  E(k)  of  (42)  has  the  same  structure  as  S(k)  except  that  for 
its  entries  the  vectors  e(k)  and  v(k)  are  used.  Of  course,  e(k)  and  v(k) 
are  only  known  statistically.  They  are  assumed  to  be  characterized  by 
the  correlations 


E le(k)e‘(J.)l  = 6 D 


(4  l.i) 


E tv(k)v^Jl)l  - 6,.D' 


(4  Ih) 


where  D and  D'  are  diagonal  matrices  of  order  m and  r respectively. 
We  shall  denote  their  diagonal  entries  as  d^^^  and  d|^  respectively. 

Further,  let  d denote  the  trace  of  D. 

Now,  the  outer  product  of  (42)  with  itself  when  summed  over 


i!L' 


k=0, 1 , . . . ,K  yields 

T T 

Y 0 Y = Y Py 


where 


>;  s‘  (k)s(k), 

k=0 

T 

Z F/(k)E(k), 
k-0 


and 


E p 


Each  of  the  identity  matrices  in  the  above  equation  is  of  order 
(N+1)  X (N+1) . Following  a procedure  similar  to  that  in  the  previous 
subsections,  it  can  be  shown  that  the  parameter  vector  y can  be  esLimated 
by  solving  the  eigenvalue  problem 

[c"^  G - e Uy  = 0 

The  eigenvector  y corresponding  to  the  minimum  eigenvalue  is  the  desired 
estimate  of  the  parameter  vector. 

2 . 4 MEASUREMENT  FILTERS 

The  problem  of  identifying  the  desired  parameter  vector,  y , was  solved 
in  the  previous  section  via  the  key  equation  (19).  The  matrices  S(k)  and 
E(k)  were  formed  by  use  of  row  vectors. 


e^(k,0)=  e^(k) 


r* 


I 


= 

= 

e(^>(k)  = 


[y^(k) , 

[Uj(k), 

[e^(k). 


, y^(k-N)] 

, Uj(k-Nj); 

, e^(k-N)] 


(49) 


( 


; ! 


f . 


li^- 

I: 


These  vectors  can  be  regarded  as  measurements  resulting  from  a cascade 
of  unit-delays  (z”^  each)  processing  the  signals  y^(k),  (k) , 

and  e^(k).  In  practice,  better  identification  results  are  achieved  if 
these  unit-delays  are  replaced  by  first  order  digital  filters  of  the  form 
q^/(l  - Qj^z“^) , where  |q^|  < 1 for  stability. 

The  precise  arrangement  is  shown  in  Figure  4-  Of  course,  the  error 
signals  e^(k),  1=0,.., m are  not  physically  available.  Rather,  they  are 
addltively  contained  in  y^fk),  i=0,..,m  respectively.  The  signals  resulting 
from  this  arrangement  are  denoted  as  y^Ck,)!,),  Uj(k,X.),  and  e^(k,2.)  ; 

£,=0,1,  .,  ,N. 

The  matrices  .S(k)  and  H:(kl  are  now  formed  by  use  of  the  row  vet  tors 

.(i) 


y (k)  = [yj(k,0) yj(k,N)] 

u^^^(k)  = (Uj(k,N-Nj) Uj(k,Nj)] 

e^^^k)  = le^(k,0) e^(k,N)] 

exactly  in  the  same  fashion  as  they  were  formed  by  use  of  the  vectors 
(35)  in  the  earlier  sections.  The  important  fact  concerning  these  new 
matrices  is  that  the  following  relation  holds; 

S(k)  A = E(k)  A 

where  the  vector  X Is  related  to  the  desired  vector  y in  a known  linear 
way.  Specifically,  it  can  be  shown  that 

y = r X 
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(50) 


(51 ) 


(52) 


■ ij 

I 


I 


i 


where  F Is  a block  diagonal  matrix;  that  Is 

r = dlag  ...  •••  ...  ...  r^} 

wherein  each  block  Fj  Is  repeated  m times.  Further,  F^  Is  a 
(Nj+l)x(Nj+l)  matrix  whose  p,  V entry  can  be  generated  by 

F = — — (F  - Q r ,) 

%-N.+N  VNj+N  P.v  -1^ 


where 


= jl  v=l 
[o  v>l  ' 


Each  of  the  formulations  in  Section  IV  can  be  restated  in  terms  of 
equation  (51)  and  the  corresponding  G and  P matrices  with  the  following 
modifications : 

• In  place  of  vector  y,  the  vector  \ must  appear  in  all  of  the 
equations  in  Section  IV. 

• The  correlation  matrix  C ( in  (31),  (37)  or  (47))  must  be 
replaced  by  the  processed  correlation 

matrix  Z. 

• After  estimating  X,  the  estimate  of  the  desired  parameter  vector 

A 

y is  computed  as  Y ” 


As  to  the  matrix  Z,  for  the  formulation  given  in  section  2.3.2  tlic  |i,v 

entry  of  Z is  given  by 

M-1  i i 

F E ( E e(i-j)  h (j)  ] ( E e(i-k)  h (k)  1 
i=0  j=0  ^ k=0  ^ 

where  h (k)  is  the  impulse  response  of  the  filter 


1 \ 


For  uncorrelated  noise  the  ii,v  entry  of  Z is 


>:  (M  - i)  h (i)  h ,(i)  . 

1=0 


I 


I 


Similar  expressions  for  the  processed  correlation  matrix  Z can  be 
derived  for  the  cases  in  subsections  2.3.3  and  2.3.4.  Indeed,  for  the 
former  case  the  expression  in  (57)  applies  except  that  h (k)  now  represents 
the  impulse  response  of  the  filter 
U-1  qp 

n r (58) 

£=1  1 - 

Nothing  has  been  said  so  far  concerning  the  choice  of  the  parameters 

and  characterizing  the  measurement  filters.  Without  loss  of 

2 

generality  let  us  choose  the  gain  parameter  as  = 1 - . As  for  the 

choice  of  we  note  that  the  measurement  filters  are  low  pass  filters  with 
cutoff  frequencies  "Ln  radian/second.  They  serve  to  decorrelate  the 

measurement  signals.  When  the  input-output  data  is  uncorrupted  by  noise, 
very  accurate  identification  is  achieved  by  a wide  range  of  choice  (over 
(0,1))  of  the  measurement  filters.  However,  when  the  data  is  corrupted  by 
noise,  as  must  be  the  case  in  real  flight  test  data,  then  the  choice  of 
becomes  the  key  to  the  success  of  accurate  identification.  Experience  has 
shown  that  a set  of  filters  with  the  cutoff  frequencies  suitably  spaced  within 
the  approximate  pass  region  of  the  system  are  highly  effective.  (For  example, 
values  that  would  attenuate  the  output  signal (s)  power  evenly  over  100  to  \0Z  — 
see  Subroutine  FINDQ) . Often  choosing  all  to  be  identical  is  adequate. 


i 


I ■ 

y 


I f 

1 I 

I i 

I 

> 

i 


[ 

r 


I 
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2.5  STATE  MODEL 

It  must  first  be  emphasized  that  even  though  we  had  Initially  formulated 
the  system  in  a state  variable  form,  the  identification  was  actually  carried 
out  in  terms  of  a transfer  function  matrix.  Let  us  denote  the  identified 
transfer  function  matrix  as  H(z)  which  when  transformed  to  an  equivalent 
s-domaln  form  becomes  H(s) . The  equivalence  may  be  impulse,  pulse  or 
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trapezoidal  equivalence.  When  the  sampling  rate  is  adequately  high  each 
of  these  equivalence  transformations  yields  approximately  the  same  matrix  H(s). 

Now,  a state  model  can  be  obtained  by  the  procedure  discussed  below. 

It  is  assumed  for  simplicity,  and  no  loss  of  generality  in  a practical 
sense,  that  the  poles  of  H(s)  are  simple.  Thus  H(s)  has  the  form 

H(s)  = (B^j(s)]  (59) 

n (s+oi.) 

^-l 


where  B^j(s)  are  the  numerator  polynomials.  A partial  fraction  expansion  yields 


n 

H(s)  = R + Z 
i=N 


s-KXj 


(60) 


The  residue  matrices  R have  the  same  dimensionality  as  H(s) . They  are 
m X r each.  Further,  these  residue  matrices  R^  can  be  decomposed  as 


i 'vi  i 
R = R L 


(61) 


'Vi 

where  R is  a matrix  whose  columns  are  simply  a maximal  set  of  linearly 
independent  columns  of  R^.  It  can  be  shown  that  if  the  system  was  identified 
perfectly,  i.e.,  if  the  identified  matrix  H(s)  had  coincided  with  the  true 
transfer  function  matrix,  then  each  of  the  matrices  R would  have  a rank  equal 
to  1.  In  practice,  this  will  not  exactly  be  the  case,  particularly  if  one  or 
more  of  the  system  modes  does  not  appear  strongly  in  any  of  t)ie  outputs.  In 
reference  [2]  such  a system  was  called  as  ill  posed. 

Let  us  assume  that  the  system  is  well-posed  so  that  eacti  of  the  matrices 
R^  lias  an  essential  rank  equal  to  1.  A matrix  will  be  said  to  have  an  esseiilial 
rank  equal  to  1 if  Its  column  (or  row)  vectors  lie  in  a two-sided  poly-cone 
of  sufficiently  small  angle.  The  analyst  must  decide  for  himself  if  the  angle 
of  the  two  sided  poly-cone  in  a particular  case  is  small  enough.  Assuming 

that  this  has  been  done  one  determines  the  decomposition  (61)  such  that  R is  a 

'Vi 

column  vector  and  the  error  is  minimized  in  some  sense.  For  example,  R 
could  be  taken  to  be  the  median  of  the  poly-cone  (see  Fig.  5). 
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1 


s+1 


0.7096 

0.7046 


[0.6789  4.0311]  + 


s+2 


s+3 


0.9530 

0.3029 

b.9910“ 

0.1338 


[-3.1054  -9.2049] 


[3.5087  6.9584] 


Note  that  the  poly-cone  angles  for  the  three  residue  matrices  were  1.6  , 
1.14°  and  0.51  . The  state  model  can  be  written  as 


dt 


’*1 

■-1 

~ 0.6789 

4.0311 

= 

-2 

^2 

+ 

-3.1054 

-9.2049 

-3 

3.5087 

6.9584 

3 

L 3. 

— 

— 

^°2^ 


'^r 

'0.7096 

0.9530 

0.9910* 

-^2- 

0.7046 

0.3029 

0.1338_ 

It  should  be  pointed  out  that  if  a pair  of  eigenvalues  turns  out  to  be 
complex  conjugate  pair,  then  the  entries  in  the  matrices  of  the  state  model 
will  appear  in  complex  conjugate  rows.  Further,  the  corresponding  state 
variables  will  also  be  complex  valued,  even  though  the  output  vector  y 
will  turn  out  to  be  real.  Such  a state  model  can  be  transformed  into  an 
equivalent  real  state  model  by  use  of  the  conversion  matrices 


J = 


'i  r 

"i  -j~ 

1 

1-1  1 

72 

.j  ~j  - 

’ ^ = 72 

_i  j_ 

Consider,  specifically,  that  the  first  two  eigenvalues  in  (62a)  arc  complex 
(and  conjugates  of  each  other).  Then  the  equations  (62a)  and  (62b),  i.e.. 


F^x  + F^u,  y 


3 4 

F X + F u can  be  transformed  into  a real  state 


model  by  defining 


P X 
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(63) 


where  P is  a matrix  whose  top  left  hand  corner  is  the  matrix  J,  the  re- 
maining diagonal  entries ’ate* unity  and  all  other  entries  are  zero.  The  real 
state  model  becomes 


X = (PF-^P~^)X  + (PF^)u 

-1  U 

y = (F  P ) X + F “ 


2.6  APPLICATION  TO  SIMULATED  FLIGHT  TEST  DATA  OF 
A SIX  MAN  SUBMERSIBLE 


The  longitudinal  dynamics  of  a six  man  submersible  was  simulated  and 
studied  with  the  objective  of  identifying  its  transfer  functions  under  non- 
zero initial  conditions.  The  transfer  function  matrix,  as  determined  by 


the  RGEORGE  computer  program  from  vehicle  geometry  and  supplied  to  us  by 


r 


the  NCSL,  was 


u(s)/6  (s) 
s 


H(s)  = w(s)/6^(s) 

0(s)/6^ (s) 


f t . / sec . /degree 
ft . / sec . /degree 


dimensionless 


0 . 001973  s^-0 . 002559s^+0 . 00003455s-0 . 00003706 
0.010383s^-0.031274s^-0. 0025350s  -0.00004588 
-0.3432s^  -0.173840s  -0.0086310 


(-0.001973) (s-0. 012012  ± j 0. 11863) (s+1 . 3208 ) 

(s  +0.0091898  ± j 0.11378) (s+0. 055798) (s+1. 4057) 

(0.010383) (s+0. 027043) (s+2.9291) 

(s  +0.0091898  ± j 0.11378) (s+1. 4057) 

(-0.3432) (s+0. 45072) 

(s+  0.0091898  ± j 0.11378) (s+1. 4057) 


where  the  common  denominator  polynomial  D(s)  is 


D(s)  = s^  + 1.47987s^  + 0.11833s^  + 0.020484s  + 0.00102196 
= (s  + 0.0091898  ± j 0.11378)(s  + 0.055798)(s  + 1.4057) 


i;; 


In  the  simulation  studies  performed,  a sampling  interval  of  0.5  sec. 
was  used.  The  equivalent  z-domaln  transfer  function  matrix,  obtained  on 
a pulse  Invai lance  basis  [6], [7],  is 

-1  -2  -3  1 

-0. 00094609+0. 0023897Z  -0.0019407z  +0.0004954z 

-0.0067869  +0.014612z"^  -0.0090902z“^+0.0012629z"^ 

-1  -2  -3 

-0.036992  +0.033999Z  +0.027069z  -0.024456z 

where 


■ WzT 


D(z)  = 1 -3.45527z“^  +4.38953z”^  -2.41136z”^  +4.77143z”‘^ 

The  response  of  the  vehicle  to  a unit  pulse  stern-plane  input,  computed 
by  use  of  H(z)  above,  is  shown  in  Figure  6.  Note  that  the  short  period 
transient  governed  by  the  mode  exp(-1.4057t)  manifests  Itself  only  in  the 
forward  velocity  u(t).  The  other  modes  manifest  in  each  of  the  three  long- 
itudinal variables  u,  w and  9.  The  problem  of  identification  is  therefore 
marginally  well-posed,  or  perhaps  ill-posed, for  fourth  order  identification. 


ise-free  study 

Three  different  runs  were  made  on  the  computer  program  MGRAM  (multi- 
variable  Gram  program  with  Initial  conditions  capability) . 

Run  1:  Non-zero  initial  conditions,  with  zero  stern  plane  input. 
Specifically, 


u(0)  = -0.018 

w(0)  = -0.022 

0(0)  = -^*31 
0(0)  = 0.0 

The  response  is  shown  in  Fig,  7.  The  system,  i.e.,  the  initial 
conditions  part  h(z)  in  (8)  (or  its  s-domain  equivalent),  was 
Identified  with  insignificant  error.  Hence,  the  details  of 
Identification  are  not  presented  (they  are  presented  for  the 
more  interesting,  noisy  data  case).  In  reconstruction,  the  ratui 


1 


:] 


I 
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of  rms  error  to  rms  signal  were 


O.OOOOU  for  u 

0.00001%  for  w 

0.00001%  for  0 

Run  2:  Zero  initial  conditions  were  assumed  and  a stern  plane  input 

0.5°  for  £ 62.5  sec 

6 (t)  » -0.5°  for  62.5  <t<  125.0  sec 

s — 

0 for  125  <t  ^250  sec 

was  applied.  The  response  is  shown  in  Fig.  8.  Again  the 
system  was  identified  almost  perfectly.  Rms  error  to  rms. 
signal  ratios  turned  out  to  be 
0.00001%  for  u 

0.00001%  for  w 

0.00001%  for  0 

The  details  are  omitted  (the  more  interesting,  noisy  data 
case  is  presented  later  in  detail) 

Run  3:  Nonzero  initial  conditions  as  in  Run  1 and  a stern  plane  input 
as  in  Run  2 were  used.  The  response  is  shown  in  Fig.  9.  Note 
that  the  response  in  this  case  is  the  sum  of  the  responses 
in  Figures  7 and  8,  since  the  system  is  linear.  As  might 
be  expected,  the  system  was  identified  with  insignificant  errors. 
The  identified  system  appeared  in  the  form 
(H^(z)  ; h(z)] 

or  the  equivalent  s-domain  form 
(H^(s)  ; h(s)l 

where  H^(s)  is  the  part  identified  in  Run  2,  while  h(s)  is  the 
initial  comlitlons  part  identified  in  Run  I. 
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Nolsy-Data  Study 

The  dynamics  of  the  six  man  submersible  is  now  studied  using  corrupted 
data.  Each  variable  was  corrupted  with  white  noise  such  that  the  ratio  of 
rms  values  of  the  noise  to  signal  was  5%.  Three  different  runs  were  made 
on  the  computer  program  MGRAM  (multi-variable  Gram  program  with  initial  con- 
ditions capability) . 

Run  4:  The  non-zero  initial  conditions  were  the  same  as  in  the  noise- 

free  Run  1.  Stern  plane  input  was  identically  zero.  To  each  of 
the  responses  an  independent  white  noise  sequence  was  added 
(rms  ratio  = 5%).  Compare  Fig.  10  with  Fig.  7.  The  results 

of  identification  are  presented  in  Table  2.1  and  in  Figure  11. 
The  latter  depicts  the  reconstructed  outputs,  comparing  them  to 
the  true  (uncorrupted)  responses.  It  should  be  emphasized  that 
the  matrix  identified  here  forms  the  second  half  of  the  matrix 
[H^(s)  ; h(s)] 

where  the  first  half  denotes  the  transfer  matrix  between  the 
stern  plane  input  and  the  longitudinal  variables.  The  errors 
in  reconstruction  are 

0.60K  for  u 

0.89%  for  w 


0.28% 


for  0 


Run  5:  Zero  initial  conditions  were  assumed  and  a stern  plane  input  as 
in  Run  2 was  applied.  The  corrupted  response  ( 5%  noise)  is 
shown  in  Fig.  12.  Compare  Fig.  12  with  Fig.  8.  The  results 
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of  identification  are  given  in  Table  2.2  and  Fig.  13.  The  latter 
shows  the  reconstructed  outputs  together  with  the  true  (uncorriipi ed) 
responses.  Note  that  the  matrix  identified  here  was  the  input- 

output  matrix  H^(s).  The  errors  in  reconstruction  are 
0.86%  for  u 

1.06%  for  w 

0.37%  for  e 

Run  6:  Nonzero  initial  conditions  as  in  Run  1 (or  Run  4)  and  a stern 
plane  input  as  in  Run  2 (or  Run  5)  were  used.  The  response 
corrupted  by  noise  (5%  on  each  variable)  is  shown  in  Fig.  14. 

The  results  of  identification  by  the  computer  program  MCRAM 
are  given  Table  2.3  and  Fig. 15.  The  latter  compares  the 
reconstructed  outputs  with  the  true  outputs.  Both  matrices 
H^(s),  representing  the  transfer  matrix,  and  h(s),  arising  from 
the  initial  conditions,  were  determined.  The  errors  in  recons- 
truction are 

0.98%  for  u 

1.31%  for  w 

0.51%  for  0 

In  each  of  the  three  noisy  runs  (Runs  4,  5 and  6)  the  wi-ak 
mode  (exp(-l .4057t) ) , present  only  the  data  of  the  variable  u,  is  missed. 

However  the  other  modes  - - poles  as  well  as  the  associated  weighting 
factors  - - are  identified  quite  accurately.  This  is  borne  out  by  the 
rms  errors  in  reconstruction,  which  was  less  than  1.5%  in  every  case.  The 
defection  of  the  micromode  could  be  Improved  by  (a)  reducing  the  sampling 
interval  and  (b)  incorporating  a slow  exponential  phase  in  the  input, 
possibly  following  tt)c  doublet. 
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TABLE  2.1 


Results  of  Longitudinal  Dynamics  Identification  of  a Six  Man  Submersible 

Non-zero  Initial  conditions:  u(0)  ■ 0.018,  w(0)  - 0.022,  0(0)  “ 0.0 

Stern  plane  Input  6 (t)  « 0 
s 

Ratio  of  Noise  to  Signal  rms  Values  ■ 5%  for  each  output  variable. 

Q Parameters:  0.890625,  0.919922,  0.943359,  0.971680 


Identified  transfer  functions 


z-domaln 


-0.0185 

+0.0539z"^  -0.0525z"^ 

-K).0170z"^ 

-0.0262 

+0.0812z"^  -0.0832z“^ 

+0.02832“^ 

-1  -2 

-3 

-4.33 

+12. 6z  -12. 2z  ^ 

+3.92Z 

-3  -•» 

-U 

-•-5.686552 

-3.69244Z  +.899341Z 

s -domain 


0.0369s^  -0.00589s^  -0.000684  -0.0000772 


0.0524s^  +0.01018^  +0. 00370s  +0.000157 


-8.668^  -1.688^ 


-0.0636s  +0.000106 


D(s)  - s^  +.212107s^  +0.02419508^  +0. 002665745s  +0.0000989037 


(s  +0.0099371  + 0.0566106  j)  (s  + 0.0544683)  ( s + 0.139284) 


Reconstruction  Errors: 


0.60%  for  u 
0.89%  for  w 
0.28%  for  0 


TABLE  2.2 


Results  of  Longitudinal  Dynamics  Identification  of  a Six  Man  Submersible 

Non-zero  Initial  conditions:  u(0)  ■ 0.0,  w(0)  • 0.0,  0(0)  » 0.0,  6 = 0.0 

Stern  plane  Input  6 (t)  ■ doublet,  duration  » 125  sec.,  amplitude  = 0.3 
s 

Ratio  of  Noise  to  Signal  rms  Values  « 5%  for  each  output  variable. 

Q Parameters:  0.933594,  0.953124,  0.970703,  0.983643 


Identified  transfer  functions 


z-domaJ.n  

-0.00125z“^  -H).00248z“^  -0.00122z“^  -0.0000231z 

-0.000585z”^  -0.0140z‘^  +0.0291z"^  -0.0145z"^ 

0.0560z“^  -0.330z“^  +0.489z“^  -0.175z”^ 

D(z)  - 1 -2.51950z"^  +1.61875z"^  -K).326244z“^  -0.425364z“^ 


s-domaln 

-0.00221s^ 

-0.00290s^ 

+0. 0000396s 

-0.0000425 

h(s) 

1 

° D(s) 

-O.OlOOs^ 

-0.0366s^ 

-0.00291s 

-0.0000533 

0.00783s^ 

-0.377s^ 

-0.200s 

-0.00982 

4 

3 . 

2 

D(s)  - s^  + 1.709628^  +0.133473s^  +0. 0236459s  +0.00116332 


- (s  +0.00917945  ±0.113694  j)  (s  + 0.0546324)  (s  + 1.63663) 

Reconstruction  Errors:  0.86Z  for  u 

1.065!  for  w 
0.37%  for  0 


h(z) 


D(z) 


degree 
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TABLE  2.3 


Results  of  Longitudinal  Dynamics  Identification  of  a Six  Man  Submersible 


Non-zero  Initial  conditions:  u(0)  = -.018,  w(0)  - -.022,  0(0)  » -4.31, 


0 ■ 0.0  Stem  plane  Input  6 (t)  “Doublet,  duratlon=125sec , 


amplitude=0.5  degree 


Ratio  of  Noise  to  Signal  rms  Values  ■ 5%  for  each  output  variable. 
Q Parameters:  0.931641,  0.953125,  0.970703,  .984131 


Identified  transfer  functions 


z-domaln  h(z) 


D(z) 


(-0.00146z“^ 

+0.00321z“^ 

-0.00202z“^ 

I (-0.0175 

+0.0283z"^ 

-0.00457z"^ 

+0.000275z"S 

1 

j 

-0.00633z"^ 

(0.000371z"^ 

-0.0159z“^ 

+0.0300z“^ 

' (-0.0179 

+0.0179z"^ 

-0.0243  z"^ 

-0.0145z“^) 

\ 

-0.0240  s'^) 

(0.0719z"^ 

-0.369z"^ 

+0.482z"^ 

I (-4.11 

+6.34z"^ 

-0.520z"^ 

-0.186Z  ) 

1 

-1.71z" 

-1 


-2 


D(z)  - 1 -2.61323Z  +1.89615z  40.0522966z  -0.335096z 

s-domaln  h(s)= 

2 


-4 


D(s) 


(-0.00271s' 


-.00376s  +0. 0000535s 
-0.0000547) 


(-0.00408s' 
3 


-0.0467s 


-0.00371s 

-0.0000682) 


( 0.0980s' 


-0.453s 


-0.257s 

-0.0126) 


(-0.0350s^  -0.0758s^  -0.0023^ 

-0.00U5) 


(-0.0348s 
3 


-0.108s 
2 


+0.04028 

+0.00238) 


(-8.21s  -18.8s 


-0.986 

+0.00929) 


D(s)  = s^  +-2.18668s^  +0.167964s^  +0.03031908  +0.00149231 


(s  +0.00926847  ± 0.113652  J)  (s  + 0.0542944)  (s  + 2.11384) 


Reconstruction  Errors: 


0.98%  for  u 
1.31%  for  w 
0.51%  for  0 
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2.7  DETERMINATION  OF  INITIAL  STATE 


Jn  the  previous  sections  It  has  been  shown  that  the  vehicle  dynamics 
can  be  Identified  successfully  In  the  form 


H(2) 


[H^(z)  ; h(z)J 


or  the  equivalent  s-domaln  form 

H(s)  = [H^(s)  ; h(s)l  (65) 

Here,  H^(s)  is  the  [nput-output  traksfer  matrix  and  h(s)  is  the  contribution 
of  the  initial  conditions.  In  fact,  if  the  continuous -time  state  model  of 
the  vehicle  were 

1 -j 

(66a) 
(66b) 


G^x  G^u(t) 


3 4 

y = G X + G u(t) 


then  analogous  to  equations  (7)  and  (8)  we  would  have 
h\s)  = G^(sl  - gS"^G^  + G* 

h(.s)  = (;’(.s1  - gS~^x(0) 


(67) 

(68) 


I' lem  Statement 

Given  H^(s)  and  h(s),  obtain  the  initial  state  for  the  state  model 
(66)  (derived  from  H^(s)  via  the  procedure  of  Section  2.5). 

Solution 

Using  the  notation  of  Section  2.5, 

N 


..1.  . " :oi  , 1 , 1 , 

H (s)  “Erl  ( — — — ) 


1=1 


s-KXj 


-b  R 


'\-i  i 0 

where  the  column  vectors  R , the  row  vectors  L , the  matrix  R and  the 

4 0 

poles  are  all  known.  Since  G = R and 

1 


.2 


L 


L L 


in  the  state  model  developed  in  Section  2.5,  we  deduce  that 


C^(8l  - G^)"^ 


, K •••♦  » •» 


s-Hx 


s-Hi 
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I 

I 

I 


Substitution  of  this  relation  Into  (68)  yields 
h(s) 


1-7^ 


s4a. 


1 'V.N 

^ R ) x(0) 


8+a 


N 

r 

1*1 


A partial  fraction  expansion  of  the  known  h(s)  and  comparison  of 


s +a. 


terms  with  the  right  hand  side  in  (69)  will  yield  the  initial  state  x(0). 


In  fact,  the  initial  state  may  be  determined  from  a single  row  of  equation  (69), 


for  example  h (s)  = h (s)  row,  provided  the  state  is  observable  in  the  variable 
3 ^ 


corresponding  to  that  row.  The  latter  is  equivalent  to  requiring  nonzero 


'i'i 


entries  in  the  vectors  R In  that  particular  row.  For  tlie  sake  of  illustration 


let  us  assume  that  sucti  Is  the  case  for  the  third  variable.  Tlien  we  liave 


N c 


3i 


n 

E 


, , ife:  ^i<»>  V 


so  that 


x^(0) 


'3i 


'V  i 

’^3 


Application 

The  importance  of  the  above  algorithm  arises  in  the  following  way.  For 
definiteness  only  the  case  of  longitudinal  dynamics  is  discussed. 

Suppose  the  vehicle  dynamics  has  been  identified  during  the  initial 
tests  as 


■ 12  3 4 

X = Gx  + Gu,  v = Gx  + Gu 


where  y = (u,  w,  0)  and  u = (6^,6  |^).  In  the  case  of  a towed  vehicle  there 
will  be  additional  inputs.  It  should  be  clear  that  this  identification  will 
require  measurement  of  each  of  the  five  variables  in  y and  u.  However, 
sucli  complete  measurements  will  be  required  one  time  only,  during  the  initial 
tests. 
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III.  ANALYSIS  OF  FLIGHT  TEST  DATA  OF  A TOWED  VEHICLE  (DSRV) 


Analysis  of  the  flight  test  data  of  the  Deep  Submergence  Research  Vehicle 
(DSRV) , supplied  by  the  Naval  Coastal  Systems  Laboratory , was  done  In 
the  following  three  stages. 

In  the  first  phase,  the  digitized  data  of  Runs  43  and  52  was  analyzed 
for  the  following  variables: 

(a)  Input  - A short  duration  (7  sec)  rudder  pulse.  It  Is  to  be 
noted  that  beyond  the  first  10  seconds  the  digitized  data  was 
clamped  to  a constant  value. 

(b)  Output  - Yaw  variable.  It  was  noted  visually  that  the  behavior 
of  this  variable  In  Run  43  was  distinctly  different  than  in  Run  52. 
Specifically,  the  response  in  Run  52  was  relatively  sluggish  sug- 
gestinig  that  Its  towing  conditions  were  different  than  for  Run  43. 

That  is,  we  feel  that  other  than  a different  rudder  Input,  there 
could  be  other  factors  causing  fundamental  differences  between  these 
runs.  Among  the  suspected  factors  are  different  stern  plane  and/or 
bow  plane  Inputs,  and  different  bridles. 

The  errors  In  identification  seemed  reasonably  low  considering  that 
the  data  were  obtained  via  digitization  of  a strip  chart  recording. 

These  errors  are  of  the  following  magnitudes. 


1 

Run 

Order  2 

Order  4 

Order  6 

43 

8.4% 

3.9% 

3.0%  i 

52 

6.8 

4.4 

3.1  ! 

r 

Graphical  display  of 

reconstructed 

outputs  can  be 

seen  in  Figures 

16  and  17.  The  identified  transfer 

functions  are 

shown  on  the 

figures. 


I 


1 i 


1 

1 

I 


In  the  second  phase  multiple-output  analysis  of  the  towed  vehicle  data 
was  carried  out  (the  approach  in  Appendix  B for  estimating  the  noise  corre- 
lation matrix  is  used).  The  variables  were 

(a)  Input  - Short  rudder  pulse  as  described  earlier 

(b)  Output  - Yaw  and  roll  variables.  The  yaw  variable  was  dis- 
cussed above.  Similar  remarks  hold  for  the  roll  variable  Inasmuch 
as  its  modes  In  Runs  43  and  52  are  distinctly  different. 

The  errors  In  identification  are  reasonably  low  considering  again 
that  the  data  were  obtained  via  digitization  of  a strip  chart 
recording. 


Run 

Order  2 (Yaw,  Roll) 

Order  4 (Yaw,  Roll) 

43 

15,8% 

9,4% 

52 

7,6% 

6,4% 

Graphical  display  of  the  reconstructed  outputs  and  the  identified 

transfer  functions,  can  be  seen  in  Figures  18  and  19. 

In  the  third  phase  analysis  of  a high  frequency  component,  clearly 
visible  in  the  Run  52  Roll,  was  carried  out.  Neither  the  fourth  nor 
the  sixth  order  identifications  had  detected  this  tremor.  Two  tech- 
niques were  used: 

Technique  1:  Identification  of  the  predomlnent  second  order  mode 
was  first  carried  out.  Call  this  component  in  roll  as  This 

was  subtracted  from  the  original  roll  data  to  yield  a residual. 

The  residual  then  used  as  the  output  data  in  a second  order 

identification. 

Technique*  2:  The  data  of  the  Run  52  Roll  variable  was  first  filUTe-d 
by  a low-pass  third  order  Chebyshev  filter,  with  a cut-off  frcquoncy  of 
0.24  hz.  Using  this  filtered  output  a second  order  identification  was 
performed.  In  order  to  reduce  the  effect  of  filter  transient,  the 
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first  four  seconds  of  the  filtered  data  were  discarded. 


I 

I 


! 

K 

f •; 

i'* 

with  both  techniques.  Identification  was  carried  out  using  two  | j 

different  inputs.  First,  the  short  rudder  pulse  presented  by  NCSL  was 
used,  since  it  could  be  presumed  that  the  tremor  was  a natural  mode 

excited  by  the  rudder  input.  Then,  a unit  pulse  located  at  t > 0 . i 

j 

was  used.  In  both  cases,  only  16  seconds  of  data  were  analyzed  in  I i 

order  to  reduce  the  effect  of  the  long  period  of  zero  input. 

Graphical  displays  of  the  results  are  found  in  Figures  20-22. 

The  following  observations  arise:  | . 

Observation  1:  Figure  20  compares  the  residual  and  filtered  output 
for  the  entire  run,  and  for  the  16  seconds  used  for  identification. 

Note  that  the  residual  behaves  in  such  the  same  way  as  does  the  ' r 

filtered  output,  both  exhibiting  a fairly  regular  pattern.  The 
conclusion  is  that  either  technique  can  reveal  the  presence  of  a 

I 

high  frequency  mode,  and  that  a tremor  characteristic  of  such  a 
mode  is  present  in  the  roll  output. 

Observation  2:  The  program  does  detect  the  tremor.  With  rudder 
input  the  poles  were  identified  as  +0.12  ± j 1.719  with  the  filtered 
output,  and  as  +0.12  t J 1.274  with  the  residual  output  (Figure  21). 

With  a unit  pulse  input,  the  poles  were  identified  as  +0.07  ± j 1.743 
with  the  filtered  output,  and  as  +0.13  ± j 1.279  with  the  residual 
output  (Figure  22).  Counting  the  number  of  peaks  in  the  tremor 
gave  an  estimated  frequency  of  about  1.8  rad/sec.  The  lower  freq- 
uency of  Che  identified  system  with  the  residual  output  is  apparently 

due  to  I he  Incl  (h.il  in  the  identification  of  the  filtered  output  !(>  j 

seconds  of  data  were  used  starting  at  time  t ■ 4 sec,  whereas  in  the 
identification  of  the  residual  output  16  seconds  of  data  starting  at 
time  c ~ 0 were  used. 


f 
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Observation  3:  Identification  of  the  filtered  output  has  very 
satisfying  aspects.  Note  the  close  relationship  of  the  reconst- 
ruction with  the  true  output,  both  In  term  of  frequency  and  phase.  ^ 

The  Initial  phase  shift  with  the  rudder  input  Is  attributed  to  the 
rudder  energy  being  concentrated  at  a tlae  somwhat  later  than  with 
the  Inpulse  input.  The  fact  that  the  technique  gave  unstable  poles  Is 
because  the  treaor  provides  a continuing,  non-decreasing  output  in 
the  interval  where  the  Input  Is  Identically  zero  (or  constant) . 

Observation  4:  The  fact  that  the  unit  pulse  input  gave  substantially 

the  same  results  as  the  rudder  input  implies  that  the  high  frequency  | 

tremor  need  not  be  a natural  mode  excited  by  the  rudder  and  could 

well  be  due  to  wave-motion.  Filtered  rudder  data  (not  8ho%m)  indicated 

the  presence  of  a tremor  of  approximtely  the  same  frequency  as  was  / 

observed  in  the  filtered  roll  data.  It  is  possible  that  this  high 

frequency  rudder  input  continued  for  the  entire  64  seconds  of  data, 

i 

but  was  not  recorded  due  to  the  limitations  of  the  digitization  I 

1 

technique  used.  This  could  have  excited  a natural  high  frequency 
roll  mode.  It  is  also  possible  that  this  rudder  motion  was  caused 
by  the  same  wave  action  that  caused  the  roll  tremor,  and  the  rudder 
actually  provided  no  further  high  frequency  input. 

I 
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IV.  UNSTABLE  VEHICLE  DYNAMICS  IDENTIFICATION 


It  is  not  uncommon  that  the  dynamics  of  the  basic  vehicle,  an  aircraft 
or  a submerged  craft,  turn  out  to  be  unstable.  Stable  behavior  must  then 
be  achieved  through  corrective  feedback.  Usually,  the  corrective  feedback 
is  provided  electronically,  but  in  the  case  of  small  diver-vehicles  it  may 
be  provided  manually.  The  situation  Is  depicted  in  Fig.  23.  The  purpose 
of  this  chapter  is  to  Investigate  the  feasibility  of  applying  the  Gram 
technique  to  the  identification  of  unstable  vehicles.  The  single-input, 
single-output  GRAM  computer  program  was  extended  to  include  the  feedback 
loop  and  the  compensator  in  digital  form.  The  program  accepts  a command 
input,  determines  the  actual  input  to  the  control  surface,  and  then  carries 
out  identification  on  the  basis  of  this  Input  and  the  vehicle  response . 
Further,  after  identification,  the  feedback  configuration  is  again  used  to 
determine  the  response  of  the  Identified  vehicle  to  the  original  command 
input  and,  in  turn,  the  reconstruction  error.  A two-man  submersible  is  used 
here  for  the  simulation  studies  performed. 


f 


Fig.  23.  Stabi I izat ion  of  Vehicle 


4.1  TWO  MAN  SUBMERSIBLE 


i The  transfer  dynamics  of  two-man  SUB,  as  predicted  by  the  RGEORGE 

program  and  supplied  to  us  by  the  Naval  Coastal  Systems  Laboratory,  is  given 
J below  for  ready  reference.  The  bow  planes  were  extended,  the  speed  was  4.5 

I 
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knots . 


Longitudinal  Dynamics 


u(8)/6g(t) 


ft. /sec. /degree 


H(s)  « w(8)/6^(t)  ; 


ft. /sec. /degree 


0(s)/6^(t) 


dimensionless 


0.0034978  -0.00039478  + 0.0002482s  -0.000018053 


D(s)  -0.04395s-’  -0.031108^  -0.0054688  -0.000067012 


-0.305758^  -0.1741928  -0.024776 


(-0. 003497) (s-K). 35432) (s-0 . 11848) (s-0 . 12297 j 
(8+0.34925) (8+0.27437) (8-0.06364) (s-0. 2327) 


(-0.04395) (8+0.41995) (8+0.013232) 
(s+0. 34925) (8-0.063/4) (8-0.23227) 


(-0.30575) (8+0.29535 

_ (s+0. 34925) (s-0.06364) (8-0.23227) 


where  the  common  denominator  polynomial  D(s)  is 


D(s)  = s^  + 0.3277138^  - 0.0739288^  - 0.019137s  + 0.0014164 


Single  variable  identification  was  performed  on  s(s)  and  0(s),  therefore 


their  transfer  functions  are  written  as  follows 


,(-0.043958^  -0.0190388  -0.000244) 

6g(s)  0^(8)' 


e(s)  = _1 (-0.30575  -0.090303) 

iS  ’(s)  D (s) 


where 


r).^(s)  = + 0. 0533428^  - 0.088564s  + 0.0051623 


•■ilL*-!.'' al  Dynamics 


“ icf-x  (0. 052868^  + 0.038108^  + 0.021775s  + 0.006405) 
6^(s)  D(s) 


* (0.05286)(8  + 0.1486  + j 0.51374)(s  + 0.42362) 

D(s;  — 


r 


i'‘ 


6^)  “ DU)  10-015728^  - 0.017968  - 0.020133] 

- (0.01572)  (8  + 0.69637)  (s  - 1.8390) 

D(8) 

where 

D(8)  - 8^  + 0.928238^  + 0.534448^  + 0.183148  + 0.023489 


X’ 


4.2  STUDY  OF  w(8)/6  (8) 

8 

The  feedback  eetup  l8  ehown  in  Fig.  24.  The  value  of  the  gain  used  to 
stabilize  the  system  was  determined  by  a computer  program.  The  closed  loop 
poles  for  this  gain,  K * -0.15,  turned  out  to  be 
|z|  « 0.97289,  0.97618,  0.99057,  0.99057 


The  results  of  identification  for  two  different  runs  are  given  in  Tables  4.1 
and  4.2.  In  Run  1 the  command  input  was  a square  wave  with  a period  of 
10  sec,  followed  by  a slow  exponential  (exp  (-0.2t)).  In  Run  2,  the  command 
Input  was  a unit  step.  In  either  case  the  response  variable  data,  w(k) , was 
nirriipled  by  VI  (rms)  noise.  It  is  clear  from  the  results  dial  the  iirsl 
iii|iul  enables  suceessiul  Identification  (see  Tin.  25)  wblU'  I bi‘  ste|»  Input 
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i 


does  not.  Several  other  coimand  Inputs  were  also  attempted  for  which  the 
results  are  not  given.  In  particular,  a doublet  with  10  sec.  duration 
yielded  successful  identification. 

4.3  STUDY  OF  0(e) /6  (s) 

s 

The  feedback  setup  is  shown  in  Fig.  26.  The  compensator  in  the  forward 
path  was  first  designed  by  root  locus  method  in  the  s-plane.  This  compensa- 
tor turned  out  to  be  s/(8+2).  The  value  of  the  gain  (in  the  feedback  path) 
was  determined  by  a computer  program.  The  closed  loop  poles  for  this  gain, 

K = -3.5,  turned  out  to  be 

\z\  = 0.91449,  0.91449,  0.99079,  0.99079 


stern  plane 
command 


6 (k 

J 

1. 

e(z) 

i-Q.a2zi’ 

Fig.  26.  Stabilization  of  pitch 

The-  results  of  Idcntlf ication  are  presented  only  for  one  run.  Run  3.  In 
this  run,  the  command  input  was  a doublet  of  duration  10  sec.  It  is  clear 
from  Table  4.3  and  Fig.  27  that  the  identification,  as  judged  from  the  poles 
and  the  reconstruction,  is  successful.  Other  runs  performed  included  a 
command  input  as  in  Run  1 (section  4.2)  and  a step  command.  The  former 
resulted  in  excellent  identification  while  the  latter,  as  might  be  expected, 
did  not. 


V’ 
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4.4  Study  of  Lateral  Dynamics 

Although  the  lateral  dynamics  Is  stable,  identification  runs  were 
performed  on  the  roll  and  yaw  rate  variables.  The  results  of  identification 
are  presented  in  Table  4.4  and  4.5.  The  rudder  input  in  each  of  these  runs. 
Run  4 and  Run  5,  was  the  same  as  in  Run  1.  That  is,  the  rudder  was  deflected 
in  a square-wave  manner,  with  10  sec.  period,  followed  by  an  exponential 
restoration  (exp  (-0.2  t)).  The  results  of  reconstruction  are  shown  in 
Figures  28  and  29. 

Discussion  - From  the  experiments  performed  it  appears  that  the  Gram 
technique  can  be  used  to  identify  the  dynamics  of  an  unstable  vehicle . To 
insure  the  success  of  the  procedure,  the  unstable  vehicle  must  be  embedded 
in  a corrective  feedback  loop  such  that  the  feedback  system  is  stable.  The 
input  and  output  of  the  vehicle  Itself  are  then  used  for  the  identification 
procedure.  This  requirement  is,  in  fact,  quite  mild  Inasmuch  as  during  the 
actual  tests  the  vehicle  is  made  stable  through  pilot  feedback,  or  through 
automatic  feedback.  The  next  point  to  note  is  the  following.  After  identi- 


fication has  been  performed,  the  reconstructed  output  must  be  computed  by 


use  of  the  same  feedback  system  as  used  before.  This  will  be  easily  possible 


escription  is 


available.  However,  when  the  feedback  is  manual,  such  description  will  be 


unavailable;  therefore,  it  wll 


ossible  to  compute  the  reconstruction 


Reconstruction  output  computed  on  an  open  loop  basis  can  be  e 


diverge  from  the  measured  output,  eve 
fled  reasonably  well. 


vehicle  parameters  are  identi- 
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TABLE  4.1 


Results  of  Vehicle  Dynealcs  Identification  for  a Two  Man  Subaerslble 


Analysis  of  w/5 


Conoand  Input 


0 ■ 0.0 
w 


o - 0.05 

V 


w(s)  _ -0.043163  s - 0.020603  s -0.00024805 

8^  + 0.0898481  8^  -0.0984324  s + 0.00578694 


Poles:  -0.384763 
0.0655826 
0.229333 


Reconstruction  Error  (RMS)  0.48Z  for  Q -0.99 


Other  values  of  Q yielded  larger  errors: 


1.966Z  for  Q - 0.96 


0.563X  for  Q - 0.98 


TABLE  4.2 


Results  of  Vehicle  Dynaalcs  Identification  for  a Two  Man  Submersible 


Analysis  of  v/6^ 
Command  Input 


99.8 


a - 0.0 

w 


- 0.05 


w(s)  _ -0.015291  s^  - 0.31351  a 0.41086 

s^  + 3.02462  s^  -2.51836  s -4.15994 


Poles:  -0.931123 
1.31193 
-3.40542 


Reconstruction  Error  (RMS) : 


7713*  for  Q - 0.96 


TABLE  4.3 


Results  of  Vehicle  Dynamics  Identification  for  a Two  Man  Submersible 


Analysis  of  9/5^ 


9(s)  _ 0.002777738  s^  -0.302842  s -8.44054 

s^  + 0.0326661  s^  -0.0838420  s + 0.00487453 


Poles  : -0.330748 
0.0625812 
0.235501 

Reconstruction  Error (RMS):  0.754Z  for  Q ■ 0.99 


Other  values  of  Q yielded  larger  errors:  0.780%  for  Q > 0.96 

0.758%  for  Q - 0.98 
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Results  of  Vehicle  Dynaalcs  Identification  for  a Two  Man  Submersible 


<f>(s)  ^ 0.0301947  3^  - 0.0324445  s^  + 0.145450  s -0.218783 

s^  + 6.04558  s^  + 3.05609  s^  + 1.67827  s + 0.246480 

Poles:  -0.154865  ± 0.460733  j 
-0.188050 
-5.54780 

Reconstruction  Error  (RMS):  2.15Z  for  Q - 0.96 


Other  values  of  Q Yielded  larger  errors:  2.970%  for  Q 0.98 

5.483%  for  Q - 0.99 


TABLE  4.5 


Results  of  Vehicle  Dynamics  Identification  for  a Two  Man  Submersible 


6 (s) 
s 


-0.0519855  s^  -0.0375124  s^  --0.0198798  s -0.00598798 
8^  + 0.924383  s^  + 0.495825  + 0.173456  s + 0.0213749 


Poles:  -0.133581  ± 0.455095  j 
-0.213805 
-0.444416 


Reconstruction  Error  (RMS):  1.49%  for  Q ■ 0.96 

I 

1 


Other  values  of  Q yielded  larger  errors:  3.97%  for  Q = 0.98 

5.07%  for  Q = 0.99 


K f/: 
h P- 
» C. 


l:  t 


<•)  Response  vs  stern  plane 
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'• 
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c 
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-6  cr 
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RfcS'’GN$f  VS  RtCCN'jT'RUC'TeD  RESPONSE: 

; I . I I . . i.lJj..-.  ij_l_l_lll.^-l-l-Ll-l-l- 

i 

1 


(b)  Reconstructed  output 

FIr.  25.  IdrntlfJratlon  of  an  unstable  vehicle: 

0 “ 0,  0 - 0.05.  A • 0.1,  Q • 0.99 


lOKliUrTLO  INPUT  a'.,i  OUTIHI 


e(.)«  <•) 


Response 


<«}  Response  vs  stern  plsne 


TfiLT  'iflS'-ON'-’F  V?  FECOUP-fRUCTFD  fiESFONGF 


e(.)/6  (.) 


True  output 


Reconstructed 


(b)  Reconstructed  output 


Fig.  27.  Identification  ot  an  unstoble  veMcle 
0 • 0.  0 • 0.05.  f.  • O.l,  0 • 0.9« 


tj'.'a'ruu  TN^'iT  hsn  ooTPur 


(•)  Response  vs  rudder 


TRUE  RESF^ONSF  VS  RFCONSTRUCTFO  RESf''0\SE 

f 1 1 I 1 1 j 1 1 . 1 t-l  i-  1-^  i J-<  < « 1 » < i aJ_I  .i,!  i iJ^L^  ,1  -Lj  ijjJj.li._l. 


C lU/ 


(b)  kuconstrurted  output 

Fig.  29>  Identification  of  an  unstable  vehicle: 

0^  - 0.05,  0^  - 0.05,  A - O.I,  Q - 0.96 
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V.  A TWO  STAGE  IDENTIFICATION  SCHEME:  GRAM-TAYLOR 


As  stated  in  Chapter  I,  the  Gram  technique  does  not  require  any  Initial 
estimates  of  the  vehicle  parameters;  additionally.  It  Is  noniterative. 

The  latter,  while  being  an  Important  advantage,  points  to  the  fact  that  the 
method  lacks  the  ability  to  Improve  results  obtained  through  a single  pass. 

On  the  other  hand,  the  NASA  program  developed  by  Taylor  and  Illff  [3} 
suffers  from  somewhat  complementary  defects.  Specifically,  Its  Iterative 
algorithm  does  not  always  converge  (to  the  true  values)  when  the  Initial 
estimates  are  not  reasonably  close  to  true  values.  It  would  therefore  appear 
that  cascading  the  Gram  and  Taylor  techniques  could  alleviate  the  short 
comings  of  either  method. 

The  cascading  of  MGRAM  to  TAYLOR  (we  shall  refer  to  the  NASA  program 
as  TAYLOR)  involved  three  things:  1)  conversion  of  the  transfer  function 
matrix  obtained  from  MGRAM  to  a state  model,  2)  Implementation  of  the 
TAYLOR  program  on  our  computer  and  3)  testing  of  GRAM-TAYLOR.  We  shall 
discuss  each  of  these  three  In  brief. 

5.1  STATE  MODEL 

Having  Identified  the  transfer  matrix  H(s)  of  the  vehicle  by  use  of 
MGRAM,  one  can  obtain  the  state  model  by  the  procedure  given  in  Chapter  11 
(Section  2.5).  An  alternative,  and  somewhat  simpler,  procedure  is  available 
for  the  case  where  the  Input  Is  a single  variable.  Specifically  consider 
that 

y(s)  “ H(s)  u(s)  (1) 

where  y Is  m-dimenslonal , us  Is  scalar  and  H Is  the  m-dlmenslonal  column 
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(a)  If  the  measured  outputs  Included  all  of  the  state  variables,  the 
program  converged  and  found  all  of  the  unknown  parameters  even 
when  the  measurements  were  corrupted  with  noise  (5Z  rms) . The 
Initial  estimates  were  all  taken  to  be  zero. 

(b)  Vrhen  the  measured  output  consisted  of  only  one  of  the  state 
variables,  specifically  the  pitch  angle,  the  algorithm  failed  to 
converge. 

5 . 3 GRAM-TAYLOR  TESTS 

The  longitudinal  dynamics  of  a Six  man  submersible,  considered  earlier 

In  chapter  2,  was  studied.  The  Initial  conditions  and  the  stern  plane  Input 

were  the  same  as  In  Run  5 of  chapter  11,  namely 

6 (t)  »,  0.5°  for  0 < t < 62.5  sec 

s - - 

<-0.5°  for  62.5  < t 1125.0  sec 
0 for  125  < t 1 250  sec 

u(0)  - w(0)  - 6(0)  - 6(0)  - 0 

The  output  measurements  u,  w and  6 were  corrupted  by  5%  rms  noise.  Identifi- 
cation via  MGRAM  yielded  the  same  results  as  In  Table  2.2  repeated  In  Table 
5.1  for  convenience.  The  state  model  was  developed  In  the  form  discussed 
on  the  previous  page  specifically,  the  matrices  of  equation  (2)  turned  out 
to  be 


.-6.0011633  -0.023646  -0.13347  -1.7096.*  1 1 _ 


] 

L 
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V 


I 


-0.0000425 

+0.0000396 

-0.00290 

-0.00221 

F » 

-a  0000533 

-a  00291 

-0.0366 

-0-0100 

-0.00982 

-0.200 

-0.377 

+0.00783 

The  vector  of  parameters  that  the  TAYLOR  program  was  given  to  Improve 
upon  was  a 20  -*  dimensional  vector.  It  consisted  of  the  four  bottom  entries 
of  A,  and  all  of  the  twelve  entries  of  the  matrix  F.  Since  none  of  the 
state  variables  is  measured  directly,  the  weighting  matrix  for  the  output 
error  in  the  algorithm  [3]  is  taken  to  be 


0 0 0 
0 0 0 
0 0 0 


0 

0 

0 


0 0 
0 0 
0 0 


D1 


0 

0 


0 

0 


0 

0 


0 

0 


0 

2173 


0 

0 


0 

0 

0 

0 

0 


0 0 0 0 


0 82 


0 


0 


0 


0 


0 


0 0.018J 


The  result  of  iterative  improvement  by  the  TAYLOR  program  are  presented 
in  Table  5.2. 

A second  experiment  was  performed  on  the  simulated  data  of  the  six  man 
submersible.  In  this  experiment  the  nonzero  initial  conditions  and  the 
stern  plane  input  were  chosen  to  be  the  same  as  in  Run  6 of  chapter  II, 
namely 


fi^(t) 


0.5° 

for 

0 

< 

t 

< 

62.5 

sec 

-0.5° 

for 

62.5 

< 

t 

< 

125 

sec 

0 

for 

125 

< 

t 

< 

250 

sec 
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and 


u(0) 

- -0.018 

w(0) 

- -0.022 

6(0) 

- -4.31 

6(0) 

= 0.0 

Identification  via  MGRAM  yielded  the  same  results  as  in  Table  2.3  of 
chapter  II,  repeated  here  in  table  5.3  for  convenience.  The  state  model 
was  developed  as  in  the  first  experiment,  specifically, 


■ 0 

1 

0 

0 

"o  " 

0 

0 

1 

0 

0 

. b = 

0 

0 

0 

1 

0 

^0.0014923 

-0.030319 

-0.16796 

-2.1867 

- 1 - 

-0.0000547 

+0.0000535 

-0.00376 

-0.00271 

-0.0000682 

-0.00371 

-0.0467 

-0.00408 

-0.0126 

-0.257 

-0.453 

40.0980 

As  in  the  first  experment,  a 20  - dimensional  vector  of  Initial  esti- 
mates was  formed  and  used  in  the  TAYLOR  program.  Of  course  the  stern  plane 
input  and  the  measured  outputs  (u,  w and  6;  each  corrupted  by  5%  rms  noise) 
was  also  supplied  to  the  program.  The  results  of  iterative  improvement  are 
shown  in  Table  5.4. 

Discussion:  In  each  of  the  two  experiments  above,  the  improvement  of 
the  apriori  parameter  estimates  Obtained  from  MGRAM)  was  not  appreciable 
although  convergence  of  the  Newton  Method  did  occur.  Lack  of  improvement 
in  the  identification  of  the  micro-model  is,  as  in  the  case  of  MGRAM, 


.1 

.1 


I j 


ft 


12 


a 

largely  due  to  the  fact  that  the  slow  sampling  rate  and  125  sec  doublet 

stern  plane  are  not  able  to  make  this  mode  observable.  Recall  this  | 

mode  is  present  only  in  the  variable  u.  Typical  computer  run  times 
on  the  IBM  360  were  250  sec  for  MGRAM,  and  800  sec.  for  TAYLOR. 
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TABLE  5.1 


Results  of  Longitudinal  Dynamics  Identification  of  a Six  Man  Submersible 


Non-zero  Initial  conditions:  u(0)  « 0.0,  w(0)  ■ 0.0,  0(0)  - 0.0,  6 = 0.0 

Stern  plane  Input  <5  (t)  « doublet,  duration  « 125  sec.,  amplitude  •»  0.5  degree 
s 


Ratio  of  Noise  to  Signal  rms  Values  • 5!S  for  each  output  variable. 
Q Parameters:  0.933594,  0.953124,  0.970703,  0.983643 


Identified  transfer  functions 


z-domaln 


h(z) 


D(z) 


-0.00125Z 


-1 


-K).00248z"^  -0.00122z"^  -0.0000231z  ^ 


-0.000585z“^  -0.0140z"^  +0.0291z"^  -0.0145z~^ 


0.0560Z 


-1 


-0.330Z 


-2 


+0.489Z 


-3 


D(z)  - 1 -2.51950z~^  +1.61875z"^  +0.326244z”^  -0.425364z 


-0.175Z 
-4 


-4 


s-domaln 


h(s) 


D(s) 


-0.00221s^  -0.00290s^  +0. 0000396s  -0.0000425 


-O.OlOOs^  -0.0366s^  -0.00291s  -0.0000533 


0.00783s^  -0.377s^ 


-0.200s 


-0.00982 


D(s)  - s^  + 1. 709628^  +0.133473s^  +0. 0236459s  +0.00116332 


(s  40.00917945  ±0.113694  j)  (s  + 0.0546324)  (s  + 1.63663) 


Reconstruction  Errors: 


0.86%  for  u 
1.06%  for  w 
0.37%  for  6 


li 
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TABLE  5.2 


Results  of  Logitudinal  Dynamic  Identification  of  Six  Man  Submersible  via 
TAYLOR  program  using  estimates  of  MGRAM  for  initialization 

9 

Zero  Initial  conditions:  u(0)  = 0,  w(0)  = 0,  6(0)  = 0,  6(0)  = 0 
Stern  Plane  Input:  6 (t)  = doublet,  duration  = 125  sec,  amplitude  =0.5° 
Ratio  of  Noise  to  Signal  rms  Values  = 5%  for  each  response  variable. 
Convergency  reached  in  10  iterations. 


Identified  State  Variable  Matrices 


A = 


0 

0 

0 


1 

0 

0 


0 

1 

0 


0 

0 

1 


-0.0011659  -0.023491  -0.13365  -1.6961J 


“0.0000425  +0.0000416-0.00291 
-0.0000533  -0.00289  -0.0362 

-0.0098385  -0.19904  -0.370 


-0.00215 

-0.00830 

+0.00141 


Characteristic  Polynomial: 


D(s)  = sVl.7s\o.l34s^+0.0235s^+0.00117 


Error : 


= (s+0. 055218) (s+1. 6224) (s+0. 0092554+ jO . 1 1370) 

2 ^ ^ 

z = / II  observed  y (k)-reconstructed  y(k)|| 

k=l  D1 

2 

Apriori  e = 0.00744 

10th  Iteration  = 0.00739 


Reconstruction  Error:  (Reconstructed  y(k)  vs  observed  y(k)) 


4.92  / 

; for 

u 

4.89  / 

1 for 

w 

5.08  / 

for 

e 
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TABLE  5.3 


Results  of  Longitudinal  Dynamics  Identification  of  a Six  Man  Submersible 

Non-a^ro  Initial  conditions:  u(0)  •»  -.018,  w(0)  “ -.022,  0(0)  - -4.31, 

0 - 0.0  Stem  plane  Input  6 (t)  “Doublet,  duratlon=125sec , 

® amplltude=0.5  degree 

Ratio  of  Noise  to  Signal  rms  Values  = 5%  for  each  output  variable. 

Q Parameters:  0.931641,  0.953125,  0.970703,  .984131 

Identified  transfer  functions 

z-domaln  h(z)  ■ 


(-0.00146z“^ 

40.00321z~^ 

-0.00202z“^ 

1 

(-0.0175 

40.0283z"^  -0.00457z"^ 

40.000275z"^) 

1 

1 

-0.00633z"^ 

(0.000371z'^ 

-0.0159z"^ 

40.0300z"^ 

1 

(-0.0179 

40.0179z“^  -0.0243  z"^ 

-0.0145z"^) 

1 

1 

-0.0240  s"^) 

(0.0719z"^ 

-0.369z"^ 

40.482z"^ 

1 

1 

(-4.11 

46.34z~^  -0.520z"^ 

-0.186Z  ) 

1 

-1.71z" 

-1 

-2  

-3  . 

-4 

D(z)  - 1 -2.61323z"^  +1.89615z“^  +0.0522966z"^  -0.335096z 


8-domaln  hCs)* 


(-0.002718^ 

-.003768^ 

40.0000535s 

-0.0000547) 

1 (-0.03508^ 

1 

1 

-0.07588^ 

-0.00238s 

-0.00115) 

D(iT 

(-0.00408s^ 

-0.0467s^ 

-0.00371s 

-0.0000682) 

' 3 

1 (-0.0348s 

1 

-0.108s^ 

40.0402s 

40.00238) 

( 0.0980s^ 

-0.453s^ 

-0.2578 

-0.0126) 

( (-8.21s^ 

1 

-18. 8s^ 

-0.986 

40.00929) 

D(s)  - 8^  +2.186688^  40.1679648^  40.03031908  40.00149231 

- (s  40.00926847  ± 0.113652  J)  (s  4 0.0542944)  (s  4 2.11384) 


Reconstruction  Errors:  0.98Z  for  u 

1.31Z  for  w 
0.51Z  for  0 
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TABLE  5.4 


Results  of  Longitudinal  Identification  of  Six  Man  Submersible  via  TAYLOR 
program  using  estimates  of  MGRAM  for  initialization. 

Non  Zero  Initial  Conditions:  u(0)=0.018,  w(0)=0.022,  e(0)=-4.31,  e(0)=0 
Stern  Plane  Input:  6 (t)  = doublet,  duration  = 125  sec,  amplitude  = 0.5° 
Ratio  of  Noise  to  Signal  rms  Values  = 5%  for  each  response  variable. 
Convergence  reached  in  6 Iterations. 


Identified  State  Variable  Matrices 


0 

1 

0 

r 

o 

0 

0 

1 

0 

A = 

0 

0 

0 

1 

_-0. 001 5010 

-0.030043 

-0.16788 

-2.1623_ 

”-0.0000548 

+0.0000614 

-0.00376 

-0.00229 

C = 

-0.0000684 

-0.00372 

-0.0460 

-0.00495 

-0.0127 

-0.255 

-0.44811 

-0.00125 

1 

Characteristic  Polynomial: 

D(s)  ^ 

= s^  + 2.16s^ 

+ 0.168s^  + 

0.030s  + 0.0015 

= (s+0. 055191) (s+2.0887) (s+0.0092310+j0. 11373) 


Error : 


f: 


2 


6th 


Apriori 

Iteration 


2 

e 

2 


€ 


K 

= ) I I observed  y(k) 
k=l 

= 0.164 
= 0.00739 


reconstructed  y(k) 


2 

1)1 


Reconstruction  Error:  (reconstructed  y(k)  vs  observed  y(k)) 

4.92%  for  u 
4.89%  for  w 
5.08%  for  0 
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VI.  LINEAR  MODEL  FOR  A MILDLY  NONLINEAR  SYSTEM 
Linear  Identification  of  the  dynamics  of  a six  man  submersible  is  studied 
while  the  maneuvers  include  the  contribution  of  nonlinear  terms.  Specif ii:al ly 
the  purpose  of  this  chapter  is  to  study  the  degradation  in  identification  of 


the  linear  model  when  the  input-output  data  actually  pertains  to  a maneuver 
in  which  the  nonlinear  force,  or  moment,  terms  make  a non-negllgible  contri- 
bution. It  is  assumed,  however,  that  the  linear  terms  are  dominant.  To 
simulate  the  vehicle  maneuvers  required  for  this  study  the  Navy  program 
DYNAMIC,  supplied  to  us  by  the  Naval  Coastal  Systems  Laboratory,  was  imple- 
mented on  the  IBM  360/65.  This  computer  program  predicts  underwater  veliicle 
responses  by  solving  the  equations  of  motion  [9].  Two  sets  of  stability 
derivatives  pertaining  to  a six  man  submersible  were  supplied  to  us  (see 
Appendix  C) : 

Set  A Linear  Coefficients  only 

Set  B Linear  and  important  nonlinear  terms.  The  latter  were  X , Z , 

c wq  ww 

Z 1 I , M I I , Y I 1 and  N i i 
w|w|  w|w|  v|v|  v|v| 

The  maneuver  chosen  was  a depth  change  maneuver  starting  from  a steady  stati- 
level  flight  at  200  feet,  and  descending  to  300  feet  via  stern  plane  actu.-ition 
only.  The  stern  plane  deflection  was  controlled  automatically  using  pitih, 
pitch  rate,  depth,  and  depth  rate  feedback.  This  capability  in  DYNAMIC, 
allowed  u,  w,  and  6 to  be  maintained  within  tolerable  limits  so  that  tlie 


differences  between  the  responses  using  Set  B 


of  coefficients,  con- 


taining the  nonlinear  stability  derivatives,  and  Set  A,  containiT\g 

only  linear  terms,  were  mild.  The  values  of  the  pertinent  parameters  used 
bv  DYNAMIC  to  execute  the  maneuver  are  given  liere. 
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I AUTO  = 

2 

ZORD  = 

300.0 

THORD  = 

0.0 

DSMAX  = 

5.0 

ZSEN  = 

-0.02 

THSEN  = 

-0.15 

DSDOT  = 

66.7 

ZDSEN  = 

0.05 

THDSEN  = 

0.6 

The  single  input,  three  output  representation  of  the  linear  system  to  be 
identified  by  MGRAM  is  described  in  the  transfer  function  equation. 


u(s) 

1 

w(s) 

= ' h2(s) 

1 

_e(s)_ 

i_  h^Cs) 

A sampling  interval  of  1.0  second  is  used  and  the  total  number  of  samples 
in  each  signal  is  181.  Two  experiments  were  performed  where  the  first 
assumed  a system  order  equal  to  4 and  the  second  a system  order  equal  to  6; 
all  else  remained  the  same  between  the  two  experiments.  In  each  of  the  two 
experiments,  two  runs  were  made:  one  corresponding  to  the  response  based  on 
the  stability  derivative  Set  A and  the  other  based  on  Set  B.  Tiius  four  runs 
were  made  as  follows: 

Run  1:  Stability  derivative  set  A,  order  of  identification  4.  The 

results  of  identification  are  shown  in  Table  6.1.  The  stern  plane 
input  for  the  desired  maneuver  is  shown  in  Fig.  30.  The  vehicle 
response  and  the  reconstructed  output  are  shown  in  Fig.  31. 

Run  2:  Stability  derivative  set  B,  order  of  identification  4.  Tlie 

results  of  identification  are  shown  in  Table  6.1.  The  stern  plane 
input  for  the  desired  maneuver  is  shown  in  Fig.  32.  The  veliiile 
response  and  the  reconstructed  output  are  sliown  in  Fig.  33. 

Run  3:  Stability  derivative  set  A,  order  of  Identification  6.  The  results 
of  identification  are  shown  in  Table  6.3.  The  stern  plane  input 
for  the  desired  maneuver  was  the  same  as  In  Run  1 , Fig.  30.  The 
vehicle  response  and  the  reconstructed  output  arc  sltown  in  Fig.  14. 
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Run  4:  Stability  derivative  set  B,  order  of  identification  6.  The  results 


i 

j 


of  identification  are  shown  in  Table  6.3.  The  stern  plane  input 
for  the  desired  maneuver  was  the  same  as  in  Run  1,  Fig.  32.  The 
vehicle  response  and  the  reconstructed  output  are  shown  in  Fig.  35. 
Discussion:  A unit  julse  response  of  the  vehicle,  obtained  from  program  DYNAMIC 
and  using  the  coefficients  of  set  A,  had  earlier  indicated  an  unstable  longi- 
tudinal response  (in  the  pitch  variable).  Each  of  the  four  runs  confirmed 
this  fact  in  that  MGRAM  Identified  this  unstable  mode.  However,  the  recon- 
struction errors  of  the  4th  order  models,  in  Runs  1 and  2,  were  undesirable 
feee  Fig.  31  and  Fig.  33).  Inspection  of  the  responses  indicates  the  presence 
of  two  natural  frequencies;  one  appearing  in  both  w and  6 and  the  second  in 
u (Fig.  31).  Hence  two  complex  pole  pairs  are  needed.  In  addition,  two 
decaying  real  exponentials  are  observed  in  the  responses  u and  0,  one  in  u 
and  the  other  in  6.  These  visual  observations  indicated  that  a 6th  order 
model  should  be  attempted.  The  improvement  in  the  6th  order  model  over  the 
4th  model  is  evident  in  Fig.  33.  Futhermore,  the  closeness  of  the  linear 
identii..ed  model  for  the  vehicle  with  coefficients  from  set  A and  the  linear 
model  identified  for  the  vehicle  with  coefficients  from  set  B indicates  that 
linear  identification  of  dominantly  linear  system  containing  nonlinear  con- 
tributions  is  indeed  viable  and  a practical  means  of  characterizing  the 
system  (compare  Tables  6.3  and  6.4). 
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TABLE  6.1 


Results  of  Longitudinal  Dynamic  Identification  of  a Six  Man 
Submersive  by  MGRAM  - u,  w,  and  6 responses  were  generated  by 
program  DYNAMIC  using  Linear  Coefficients  only. 

Identification  Order  N = 4 

Depth  Change  Maneuver  from  200  to  300  feet. 


Identified  transfer  functions 


h(s) 


D(s) 


-O.OOllOs^ 

-0.000113s^ 

-0.0000251s 

-0.00000142 

1 

-0.0266s^ 

-0.00170s^ 

-0.000261s 

-0.00000362 

D(s) 

-0.268s^ 

-0.153s^ 

-0.0105s 

-0.00126 

+ 0.025793s^  + 0. 

020030s^  + 0. 

00049584s  + 0. 

0001 2 108 

(s  + 0.031575  ± j 0.090207)(s  - 0.018679  ± j 0.11361) 


Reconstruction  Error: 


34.7%  for  u 


8.5%  for  w 
3.3%  for  6 
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TABLE  6.2 


Results  of  Longitudinal  Dynamic  Identification  of  a Six  Man 
Submersive  by  MGRAM  - u,w,  and  6 responses  were  generated  by 
program  DYNAMIC  using  Linear  and  Nonlinear  Coefficents. 

Identification  Order  N = 4 

Depth  Change  Maneuver  from  200  to  300  feet. 


Identification  transfer  functions 


-0. 

1 


-0.00143s^  -0.000246s^  -0.0000223s  -0.00000355 


-0.0271s^  -0.00120s^ 


0.176s- 


-0.104s 


-0.000316s  -0.00000619 


-0.00420s  -0.00113 


D(s)  = s^  -0.0021104s^  + 0.019597s^  - 0.000074644s  + 0.00010438 
= (s  + 0.014577  ± j 0.10449) (s  - 0.015632  ± j 0.095565) 


Reconstruction  Error: 


47.4%  for  u 


10.2%  for  w 


8.9%  for  0 


TABLE  6.3 


I 

t 

( 


I 

i! 


I 


I 


i: 


i: 


Results  of  Longitudinal  Dynamic  Identification  of  a Six  Man 
Submersive  by  MGRAM  - u,w,  and  9 responses  were  generated  by 
program  DYNAMIC  using  Linear  Coefficients  only. 

Identification  Order  N ■ 6 

Depth  Change  Maneuver  from  200  to  300  feet. 


Identified  transfer  functions 


~0.000371s^ 

-0.000234s^ 

-0.000285s^ 

-0.00000799s 

-0.0000208s^ 

-0.000000152 

* D(s) 

-0.0297s^ 

-0. 007908^ 

-0.00449s^ 

-0.0000403s 

-0.000735s^ 

-0.000000487 

-0.281s^ 

-0.212s^ 

o 3 

-0.080  8s 
-0.00405s 

-0.  0282s^ 
-0.00015] 

D(s) 


6 

s 


+ 0.25405s^  + 0.15750s^  + 0.025891s^  + 

+ 0.000327538+ 


0.0021218s^ 

0.000014630 


- (s  + 0.040930  ± j 0.35744) (s-0. 018657  ± j 0.11361) 
(s  + 0.055291) ( s + 0.15422) 


Reconstruction  Error: 


4.0%  for  u 


5.4%  for  w 
0.91%  for  U 


I 


1 


\ 
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TABLE  6.4 


Results  of  Longitudinal  Dynamic  Identification  of  a Six  Man 
Submersive  by  MGRAM  - u,  w,  and  0 responses  were  generated  by 
program  DYNAMIC  using  Linear  and  Nonlinear  Coefficients. 

Identification  Order  N ■ 6 

Depth  Change  Maneuver  from  200  to  300  feet. 


Identified  transfer  functions 


I 

i 

I 


j 


Reconstruction  Error:  2.5%  for  u 

4.9%  for  w 
1.6%  for  0 


h(s) 


D(s) 


■-0.000875s^  -0.000514s^  -0.000556s^  -0.0000303.s^ 

-0.0000112s  -0.000000223 


-0.0313s 

-0.281s^ 


-0.00998s  -0.00350s  -0.000682s 

-0.0000318s  -0.000000341 


-0.212s 


-0.0808s-' 

-0.00405s 


-0.0282a^ 

-0.000151 


D(s)  - s^  + 0.29811s^  + 0.11050s^  + 0.0231248^  + 0.00099258s^ 

+ 0.00019881S+  0.0000064645 

=■  (s  + 0.031807  ± j 0.29261) (s  - 0.015755  ± j 0.095482) 
(s  + 0.034405) (s  + 0.23160) 
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/DYNOMIC.  LINffifv  Tftt/  DfP'H  f.lKiNOi  H'ltJlllYK 

i-l  I_1  I J 1 I I 1 1 1 1.  1 1 .]  1 J l-l.i-  X J ill  1 .1  I.  I 


TJtC  IN  SFCOr^D 


Vehicle  •l•uUlt■d  via 


30.  St«rn  pl«n«  deflection 
coefficient  tet  A. 


. • IN! 'iK  fFSl/  i.h'iHl' 


true  rcwrenst 
r«con*tructed  r«tponM 


./[■VNfiHIC.  LINF-OR  TEST/  OtPTh  ChnNGF  MPf: 


true  response 
cecooetrucied  response 


— r-nCH  - /D'’''lfiMIC  LlNt:f)R  TFST/  DFPTH  CHnNC-t  f/'-'itUVcR 


true  response 
teeoAscrucced  reeponee 


l.^ns ^ ( ud I n.i ] v..  rr%}'>Ki*.>  rt-ton  trtittrd  (rmii 

order  fd«'nt  ( f U .«l  ii*n . Vidtlt  te  AfftitliieJ  vTi  irK-ff2ri«ni 


r/  OLPTM  CHfUaif  MONFLIVi  !• 


;UB  --  STN  PLN  — /DYNAMIC  NON-LlNEPiR  M 


Fig.  32.  Stern  plane  deflection.  Vehicle  einuleted  via 
coefficient  aet  B. 


t»  - ‘.'JG  --  U - ^FST'  f'[*'’H  fHONOr 

C C C“  ] i .1  i .i  j.  l.l  1 1.  1.  I L 1 1.1  I J.  J .1.-1. 1 l-.l.L.LL_L.tU-Ll  I,.'  t 1 i.  1 1 l.J  1 J 111  t 1 i I 


u(t) 

- I true  response 

reconstructed  response 


I [ \ rr-T~|T  »"rT-pri~T  i | i i i i | i"f  i t j rTT~rjT  i i f | i i i i jT-rr-i-y 
J8  C >‘•.0  C 72.0  3C.C  JC5  IM.  n:  ]50. 


6 - K^N  SU9  — w — /DVNftrirC.  LINFPR  TEST/  DEPni  CHANGE  MPNEL'^ER 


( 20: 


Lxj.j-1,-- 


V(t) 

——  true  response 

reconstructed  response 


lie  r>.0  64.0  7?  c 90  0 iCfl  l?0  IM.  If?  I50 


6 - rm  SUB  — PlfCH  — /OTNRMIC.  LINEPR  TfST/  DEPTH  CHRfJOf  MHNtU/ER 


i J- 


■ true  response 
— — reconstructed  response 


le  0 j 1 t t f J r ; f i ; 1 I » : I 1 r t t I r f f •-  r T ! I -»  I -T  1 • r J f I f r j r I T 1 I J » : I • 
t if  0 5.  M )■  ',  ...  i,'*>  :*y.  n' 

ir. 


Fie.  14.  I onr.Ilu'llti.il  f»'sp«‘n>n*  vh  rrcoivu  mri  id  fr«*r»  nh 

ordt  f IH.'UI  1 1 Irai  Ion  . V<)i|«te  vl.i  io,iM<l«ni 


SUB  — W --  /DrWMIC.  NOfi-LINEflR  TEST/  DEPTH  CHmJC-r  MMUVER 


E:=iN  Sv.5  — ‘=;TCH  '■■  /DTNfiMIC  NON-i.  INEPiR  r.EP'H  fvi’-ltl 


FIr.  JK  LonRlliidlnrtl  rpspon^c  rp^lMinMp  Tpionst  rut  ti-d  frti»  fiih 
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APPENDIX  A 


Till-  NASA  it)nipiitc*r  program  written  by  L.  W.  Taylor  and  K.  W.  11  iff  | Jj 
nsi's  a modified  Newton  Raphson  method  for  system  Identification.  This 
program  has  been  adapted  for  implementation  on  the  IBM  360/65  computer  of 
the  University  of  South  Florida.  Because  of  some  apparent  discrepencies 
between  the  thoery  and  the  program  copy  supplied  to  us  by  the  Naval  Coastal 


Systems  Laboratory,  we  have  incorporated  some  modifications.  The 
computer  program,  now  operational,  is  used  for  the  GRAM-TAYLOR  interface 


discussed  in  Chapter  V. 

To  explain  the  modifications  done  we  consider  the  computation  of  the 


gradient  of  the  output  with  respect  to  the  vector  of  unknown  parameters. 


Till'  notation  used  is  largely  the  same  as  in  (31.  Let: 


X = Ax  + Bu 


state  equation 


X + b 


Fx  + Gu  + d 


output 


x(0)  = x(0)  + q 


initial  state 


x(0)  given 

Let  c denote  the  vector  of  all  unknown  parameters.  To  use  the 

recursion  formula  to  Improve  an  estimate  of  c,  we  need  the  gradient 

V y.  For  this  purpose  let  us  differentiate  (1)  with  respect  to  c: 
c 

• 

9x  . J)x  . 3A  . 9B  ... 

+ -^x  + ^u  (Aj 

dC  dc  dc  dc 

Letting  the  subscript  k denote  the  sampled  values  of  the  variables 
(for  example  y^  • y(kA),  where  A is  the  sampling  Interval),  and  assuming 
a stair  case  approximation  for  the  continuous  time  variables  we  have  the 


following  solution  for  equation  (4): 
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On  the  other  hand,  by  direct  differentiation  of  (2)  we  have 


Simplifications  - We  now  recognize  that  the  vector  c has  the  form 


where  p denotes  the  subset  denoting  the  unknown  entries  from  the  matrices 
A,  B,  F and  G.  The  vectors  q,  b,  d were  defined  earlier  in  equations 


we  shall  compute  each  entry  on  the  right  hand  side  of  (8)  in  as  simple 


a form  as  possible 


First  Row  of  (8)  - From  (5) 


T 


^ ’ l|?  -k  * It  “k'  • 

'IV'o  ■ “ 


/^y  ^ - vA  + ^ X + ^ u 

'ap  -'k+i  ~ *^'r)p^k+i  3p  *k+i  ap  k+i 


r,,3y^,.  , 3F  3^ 

’^^ap  \+l  ■*■  3p  \+l  ’’’  'f)p  “k+1 


Iniiial  value  for  this  can  be  determined  using  (9d)  and  other 
quant  it ies. 

Second  Row  of  (8)  - From  (5)  and  (7), 


/ A / 9x\ 

^9q\+l  " ^ ^aq^k’ 


1 


(f^)  = I identity  matrix 

o<l  0 


* (h.\  (ii\  = I 

i-i-  ^ ■'i,'  '9q  ■'0 


1 
I 

^3q  'k+1  " '3q  "k 

= F (lij 

^3q  'k+l  '3q''k+l 
Third  Row  of  (8)  - From  (5)  and  (7) 


= <1,  (ili) 

\)b^k+l  ^3b^k’ 

''3b  'k+1 


(ill)  s 0 
'3b‘'o 


^Ib  ^k+1  = ° 


Fourth  Row  of  (8)  - From  (5)  and  (7), 


/dXy  jt.  / 9X\ 

^3d^k+l  “ ^ ^3d^k’ 


<lf>0  = 0 


1 


(^  ) - 0 
^3d  ^k+1 


(^  ) = 1 
^3d  ^k+1 
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(9c) 

(9d) 

(9e) 

known 

(10a) 

(10b) 

(10c) 

(lla) 

(llb) 

(llc) 

(12a) 

(12b) 

(12c) 


1 


U 

Kk.  r.  • 


OjLscussion  - Equation  (9)  through  (12)  should  have  been  used  to  compute  the 
gradient  recursively.  For  the  most  part,  the  copy  of  the  NASA 

program  supplied  to  us  did  accomplish  this  recursion  correctly.  However, 
the  following  errors  were  found  and  corrected  by  us.  They  are  discussed 
here  for  the  benefit  of  the  many  other  users  of  the  Taylor  program 


Earlier 


I .  The  program  computed 


y k+i  “ •'^+1  ^“k+i 


2.  To  compute 

earlier  program  was  using 
3 1 

(4^).  in  equation  (9e)  . 

3p  "k  2 

3.  To  compute 

3y^ 

program  used 
uation  (10c). 

4.  x(0)  = X (0)  + q - b 

;»y* 

"7.  updating  of  was  in  error 


f ied: 


=i  - 


and  the  error  function,  J,  is 
given  by 

■'i  j 


It  now  computes 


>'  k+i  ■ "k+i  * 


y k+i  ■ ■'"k+i  ‘'“k+i  * ‘‘ 

3 1 

It  now  substitutes 

into  the  right  hand  side  of  (9e) . 


3 ^ 

It  now  substitutes 

the  right  hand  side  of  (10c) 


This  is  changed  to  x(0)  = x (0)  + q 
This  is  corrected  to  conform  to  (lOh) 


ft.  Updating  of  parameters  was  modi-  Modification: 


‘^i+l  “ 

where  at  each  iteration  a is  choosen 


to  be  the  largest  number  In  the  set, 

{1,  1/2,  ...,  1/1024),  which  decreases 
the  error  function  from  the  previous 
iteration.  If  an  a cannot  be  found  in 
the  above  set,  then  convergence  is 
assuimued . Thus  at  each  Iterat  ion 
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APPENDIX  B 


AN  APPROACH  TO  UNKNOWN  NOISE  IN  MEASUREMENTS 

In  Chapter  II,  vehicle  dynamics  Identification  was  carried  out  with 
noisy  measurements.  However,  the  characteristics  of  the  noise,  or  error, 
process  were  assumed  to  be  known.  The  purpose  of  this  section  is  to  propose 
an  approach  for  the  case  when  information  about  the  noise  process  (the  vector 
e(k)  in  (Ic) , or  e(2)  in  (11))  is  lacking.  Specifically,  a method  is 
proposed  to  iteratively  estimate  the  correlation  matrix  P of  equation  (25) . 
There  is  as  yet  no  theoretical  basis  for  the  procedure,  nor  does  it  always 
satisfactorily  converge  to  a P matrix.  Indeed,  for  an  ill-posed  system 
identification  problem,  the  procedure  results  in  oscillatory  behavior.  Dra- 
matic improvements  have,  however,  resulted  in  cases  where  the  identification 
problem  is  well  posed  (l.e.  each  mode  of  the  system  appears  reasonably 
strongly  in  at  least  one  of  the  observed  outputs) . 

To  start  the  method,  the  white  noise  assumption,  as  discussed  in  sub- 
section 2.3.2  of  Chapter  II,  is  used  to  estimate  the  initial  (zeroth  iterate) 
parameter  vector  This  vector  is  next  used  to  determine  the  system 

response,  and  in  turn  the  error 

e^^^(k)  * y(k)  - y^®^(k) 

This  error  sequence  is  used  to  compute  P^^^.  The  formulas  necessary  for  the 

fkl 

purpose  are  (18).  which  define  the  matrix  E‘  and  (27).  This  estimate  of 
error  correlation  matrix  is  used  to  determine  a new  estimate  Specif- 

ically, it  is  chosen  as  the  eigenvector,  with  its  first  entry  normalized  to 

1,  corresponding  to  the  minimum  eigenvalue  solution  of 
^(IIT  ^,^11]  . ^[1]T  p[l]  ^[1] 

f 2] 

This  solution  can  be  used  to  generate  a new  error  sequence  e^  and  the 
procedure  can  be  Iterated  until  satisfactory  convergence  occurs. 
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